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Executive  Summary 


Accurate  prediction  of  turbulent  flows  at  high  Reynolds  numbers  and  with  substantial  effects  of 
separation  constitute  a  regime  of  significant  engineering  interest.  Solutions  of  the  Reynolds- 
averaged  Navier-Stokes  (RANS)  equations  are  unable  to  represent  the  complex  physics  of  most 
separated  flows  to  sufficient  accuracy.  Large-Eddy  Simulation  (LES)  provides  a  more  realistic 
treatment  of  the  separated  regions  of  a  turbulent  flow  but  is  prohibitively  expensive  when  applied 
to  whole  domains  at  high  Reynolds  numbers.  Detached-Eddy  Simulation  is  a  hybrid  method,  com¬ 
bining  RANS  and  LES,  and  attempts  to  take  advantage  of  both  techniques  in  regions  where  each 
is  accurate  and  computationally  feasible.  In  natural  applications  of  the  method,  attached  boundary 
layers  arc  entrusted  to  the  RANS  model  with  detached  regions  of  the  flow  predicted  using  LES. 

While  DES  is  now  relatively  well  understood  in  massively  separated  flows  characterized  by  thin 
boundary  layers  prior  to  separation,  an  incorrect  behavior  can  be  encountered  in  flows  with  thick 
boundary  layers  and/or  shallow  separations.  This  motivated  a  detailed  study  of  the  flow  physics 
and  grid-sensitivity  of  the  solutions  in  the  flow  over  an  Aerospatiale  A-airfoil  at  an  angle-of-attack 
of  a  —  13.3°  and  at  a  chord-based  Reynolds  number  of  2.1  x  106.  This  flow  has  been  measured  in 
separate  experiments  (29),  (30)  and  was  the  subject  of  a  coordinated  set  of  investigations  through 
the  LESFOIL  project  (3 1 ).  The  flow  is  a  useful  test  case  since  problems  can  arise  when  the  grid 
spacing  parallel  to  the  wall  (A||)  becomes  less  than  the  boundary  layer  thickness  (fi),  which  is 
possible  either  through  grid  refinement  and/or  boundary  layer  thickening.  In  this  regime,  the  grid 
spacing  is  then  fine  enough  for  the  DES  length  scale  to  take  its  LES  value  within  the  boundary 
layer.  This  lowers  eddy  viscosity  levels  below  that  which  would  be  given  by  the  RANS  model. 
Resolved  stresses  that  are  derived  from  velocity  fluctuations  (“LES  content’’)  may  not  yet  have 
been  generated  or  might  be  weak.  Additionally,  LES  content  may  be  lacking  because  the  grid 
or  time  step  are  not  fine  enough  to  fully  support  it,  and/or  because  of  delays  in  its  generation  by 
flow  instabilities.  The  lower  total  stress  reduces  the  skin  friction,  which  can  lead  to  premature 
separation  in  extreme  cases. 

To  address  these  deficiencies,  a  new  version  of  the  technique  known  as  Delayed- Detached 
Eddy  Simulation  (23)  (DDES)  has  been  developed.  In  natural  applications  of  the  method,  DDES 
over-rides  the  DES  limiter  and  maintains  RANS  behavior  in  the  boundary  layers  regardless  of  the 
grid  spacing.  The  new  version  has  been  tested  on  the  flow  over  the  Aerospatiale  A-airfoil.  To 
ensure  that  the  new  version  does  not  degrade  predictions  in  massively  separated  flows,  the  flow 
over  a  circular  cylinder  at  Reynolds  numbers  corresponding  to  1.4  x  IQ5  and  8  x  1G6  has  also  been 
analyzed,  RANS  predictions  of  the  flow  over  an  Aerospatiale  A-airfoil  at  a  chord-based  Reynolds 
number  of  2.1  x  10fi  at  two  lower  angles  of  attack  have  also  been  included.  The  DDES  results  of 
the  flow  over  the  Aerospatiale  A-airfoil  are  essentially  grid  independent  and  match  well  with  the 
RANS  results.  In  contrast,  the  DES  results  exhibit  premature  separation  occuring  in  two  of  the 


i 


grids.  The  eddy  viscosity  profiles  showed  the  depletion  of  modeled  stresses  with  the  use  of  finer 
grids  in  DES,  as  a  result  of  the  RANS-LES  interface  moving  well  within  the  boundary  layer.  The 
flow  over  the  circular  cylinder  at  a  Reynolds  number  of  8  x  106  showed  that  the  DES  and  DDES 
predictions  are  similar  with  good  agreement  between  the  statistics  using  the  two  models. 

While  the  new  method  addresses  natural  DES  applications,  the  other  primary  thrust  of  the  cur¬ 
rent  work  reported  in  this  document  focuses  on  the  transition  from  RANS  to  LES  behavior  in  DES 
along  the  main  flow  direction.  This  can  result  in  applications  because  of  streamwise  grid  refine¬ 
ment  as  required  for  resolution  of  tlow  features  or  changes  in  geometry  that  demand  streamwise 
grid  refinement.  This  part  of  the  work  is  further  motivated  by  flows  characterized  by  shallow  sepa¬ 
rations  where  RANS  models  may  not  provide  sufficient  accuracy  and  an  LES  treatment  is  desired 
in  order  to  improve  accuracy. 

The  transition  from  RANS  to  LES  has  typically  been  studied  with  evolution  normal  to  the 
wall,  as  in  the  wall-layer  modeling  work  of  Nikitin  et  al.  (72).  Whatever  the  mechanism,  as  the 
flow  enters  the  LES  region,  a  DES  solution  is  characterized  by  a  reduction  in  the  modeled  length 
scale  which  draws  down  the  eddy  viscosity,  thereby  lowering  the  modeled  stress.  This  reduction 
in  modeled  stress  requires  a  corresponding  increase  in  the  resolved  stress  near  the  transition  layer 
in  order  to  avoid  degrading  predictions  of  the  boundary  layer.  One  of  the  techniques  that  can 
be  applied  to  quickly  increase  resolved  stresses  is  to  seed  the  near- wall  flow  w'ilh  fluctuations  in 
order  to  quickly  generate  three-dimensional  structure  A  scheme  has  been  developed  and  tested 
that  generates  velocity  fluctuations  using  a  stochastic  force  added  to  the  momentum  equations  and 
subsequent  tuning  of  the  amplifications  based  on  a  prescribed  shear  stress  profile.  An  adaptive 
control  scheme  has  been  used  to  shape  the  velocity  fluctuations  to  mimic  realistic  flow  features 
and  optimized  to  achieve  realistic  Reynolds  shear  stresses  with  minimal  computational  effort. 

Before  arriving  at  the  adaptive  control  scheme,  the  control  process  was  thoroughly  investigated 
and  an  analysis  of  the  behavior  of  the  errors  using  Proportional-Integral  (PI)  controllers  and  in¬ 
tegral  controllers  was  carried  out.  These  simulations  also  revealed  that  specification  of  controller 
gains  that  ensure  stability  of  the  system  is  not  simple.  Therefore,  the  main  advantage  of  using  the 
adaptive  control  scheme  is  that  it  obviates  the  need  to  specify  these  controller  gains  that  ensure  nu¬ 
merical  stability.  The  adaptive  control  method  is  based  on  the  MIT  rule  which  is  an  adaptive  control 
methodology  proposed  by  Whitaker  (75).  Using  the  method,  the  system  is  initially  in  the  open- 
loop  configuration  with  no  control  forces  being  applied.  The  gain  values  are  then  continuously 
altered  during  the  simulation  in  an  adaptive  manner.  The  alteration  of  gains  is  based  on  a  law  that 
involves  approximating  a  gradient-descent  procedure  seeking  the  minimum  of  an  integral-squared 
performance  criterion.  The  scheme  has  been  tested  using  computations  of  turbulent  channel  flow 
at  Reynolds  numbers  based  on  friction  velocity  and  channel  halfwidth  of  5000.  Simulation  results 
show  that  the  method  proposed  is  effective  in  achieving  the  desired  results  with  the  resolved  shear 
stress  reaching  the  target  levels  at  the  control  planes. 


Simulations  were  conducted  using  coarse  and  fine  grids,  with  wall-parallel  grid  spacings  of 
approximately  I / 1 0  and  1/20  of  the  boundary  layer  thickness,  respectively.  These  simulations 
revealed  that  the  finer  grid  is  able  to  sustain  turbulence  downstream  of  the  control  planes  and 
reduces  the  error  in  the  wall  shear  stress  from  20%  to  less  than  10%  for  a  channel  with  a  stream  wise 
extent  of  4tt  and  with  the  controllers  taking  up  approximately  30%  of  the  domain  extent  in  the 
sireamwise  direction. 
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1  Background 


Predicting  turbulent  flows  constitutes  one  of  the  principal  challenges  of  physics.  The  non-linearity 
and  chaotic  nature  of  the  solutions  and  influence  of  complicating  features  such  as  the  Reynolds 
number  and  presence  of  a  solid  surface  substantially  challenge  simulation  strategies.  As  described 
by  Spalart,  there  are  two  main  challenges  in  the  prediction  of  turbulence  (49): 

•  Prediction  of  growth  and  separation  of  the  boundary  layer; 

•  Prediction  of  momentum  transfer  after  separation. 

The  most  commonly  used  simulations  strategies  range  from  solution  of  the  steady  Reynolds- 
averaged  Navier-Stokes  (RANS)  equations  to  Large  Eddy  Simulation  (LES)  to  Direct  Numerical 
Simulation  (DNS).  In  the  recent  past,  a  number  of  hybrid  strategies  have  emerged,  motivated  by 
both  the  strengths  and  limitations  RANS  and  LES.  The  most  popular  of  these  hybrid  methods 
include  Detached  Eddy  Simulations  (DES),  Very-Large  Eddy  Simulation  (VLES),  and  Unsteady 
Reynolds  Averaged  Navier-Stokes  (URANS)  solutions. 

With  the  obvious  exception  of  DNS,  simulation  strategies  for  applications  require  empirical 
input.  As  described  in  Spalart  (49),  it  appears  that  complex  modeling  strategies  sometime  advo¬ 
cated  for  RANS  approaches  offer  relatively  little  advantage  in  (be  prediction  of  the  growth  and 
separation  of  a  boundary  layer,  as  most  simple  RANS  models  have  been  calibrated  in  thin  shear 
layers.  Predicting  momentum  transfer  after  separation  is  more  challenging  and  hybrid  strategies  as 
well  as  complex  RANS  models  have  distinct  advantages  over  simple  schemes.  White  some  strate¬ 
gies  involve  extensive  use  of  empirical  inputs  into  models,  others  such  as  LES  involve  much  less 
empiricism.  In  some  strategies,  grid  refinement  merely  improves  the  numerical  accuracy  of  the 
solution,  and  in  others,  we  could  obtain  richer  turbulence  physics  by  employing  a  finer  mesh.  The 
final  selection  of  a  turbulence  strategy  to  predict  a  flowfield  is  inevitably  a  compromise  between 
accuracy  and  computational  efficiency. 

1.1  Direct  Numerical  Simulations 

In  this  approach,  which  is  free  from  any  empiricism,  the  entire  range  of  turbulent  scales  are  re¬ 
solved.  In  DNS,  the  Navier-Stokes  equations  are  solved  without  any  modeling.  DNS  is  the  most 
fundamental  technique  as  it  provides  detailed  information  about  turbulence,  which  can  be  difficult 
to  obtain  in  experiments.  Since  all  the  scales  of  turbulence  are  to  be  resolved  using  DNS,  the  grid- 
ding  requirements  are  stringent,  and  it  must  be  ensured  that  DNS  grids  are  fine  enough  to  resolve 
the  smallest  scales  while  the  domain  should  be  large  enough  to  contain  an  adequate  sample  of  the 
largest  eddies, 

In  most  turbulent  flowfields,  energy  is  supplied  via  interactions  between  the  turbulence  and 
mean  flow  gradients  to  the  largest  eddies  (that  have  length  scales  L  of  the  order  of  the  domain). 


and  in  turn  transfer  kinetic  energy  to  the  smaller  scales  until  that  energy  is  dissipated  into  heal 
at  the  smallest  scales  (the  Kolmogorov  length  scale  ■//).  The  large  eddies  are  influenced  by  the 
domain  and  are  anisotropic.  This  anisotropy  is  lost  during  the  process  of  the  energy  cascade  and 
the  smallest  eddies  arc  assumed  isotropic. 

If  the  rate  of  dissipation  of  turbulent  kinetic  energy  (A:)  is  c,  then  it  would  be  possible  to  define 
turbulent  lengthscales  and  turbulent  timescales  on  this  basis.  By  noticing  that  u  =  \/k  is  a  velocity 
scale,  the  large  eddy  length  scale  can  be  defined  as: 

_  fc3/2 

L  =  rVk  =  —  CD 

e 

where  T  is  a  timescale  for  energy  dissipation  (T  =  k/e),  which  is  the  ratio  of  the  turbulent  kinetic 
energy  to  (fc)  the  rate  of  dissipation  (t)  of  turbulent  kinetic  energy.  This  quantity  is  also  referred 
to  as  the  integral  timescale.  Since  the  smallest  scales  of  motion  have  little  energy,  it  would  not 
be  logical  to  use  k  to  define  a  timescale  for  these  scales.  As  the  smallest  scales  are  diffused  by 
viscosity,  it  is  appropriate  to  form  a  timescale  using  v  (viscosity)  and  (. 


Tv  is  referred  to  as  the  dissipative  timescale.  Therefore,  the  ratio  of  the  timescales  of  the  large  and 
small  eddies  is, 


Similar  to  the  definition  of  the  large  eddy  lengthscale,  the  lengthscales  of  the  smallest  scales  of 
motion  (r/)  can  be  defined  using  u  and  e: 


T)  = 


1/4 


(4) 


The  ratio  of  the  largest  lengthscale  (L)  to  the  smallest  lengthscale  (r/)  is  given  by: 


k3/2 

=  e3/4T;3/4 


(5) 


L/rj  =  Re3/* 


(6) 


This  implies  that  a  three-dimensional  DNS  requires  the  number  of  mesh  points  proportional  to 
(/?eV4)  i  (24).  Similarly,  the  ratio  of  the  integral  timescale  ( T )  to  the  dissipative  timescale  (T,,) 
is  a  function  of  the  Reynolds  number  (Re/2)  (24).  It  is  also  important  that  the  integration  of  the 
solution  in  time  must  be  done  with  a  time-step  (Af).  small  enough  that  a  fluid  particle  moves  only 
a  fraction  of  the  mesh  spacing  (Ax)  in  each  step.  This  restriction  is  in  the  form  of  a  Courant 
number  (C). 


_  uAf 
C  =  ——  <  1 
Ax 


(7) 
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The  total  simulation  time  is  generally  proportional  to  the  turbulence  timescale  given  by: 

L 

T  —  — 

U 

Combining  the  above  relations  and  because  Ax  should  be  of  the  order  of  77,  the  number  of  lime- 
integration  steps  must  be  proportional  to  LjCr).  Since  Ljrj  =  Re)) one  can  estimate  that  the 
number  of  floating  point  operations  required  to  complete  the  simulation  is  proportional  to  the 
mesh  points  and  the  number  of  time-steps.  Thus,  the  CPU  time  grows  as  Re)  (25).  Due  to  the 
prohibitive  computational  cost,  the  application  ofDNS  to  complex  flows  at  high  Reynolds  numbers 
is  problematic.  However,  DNS  remains  a  useful  and  powerful  technique  to  compute  low  Reynolds 
number  flows,  and  provides  pathways  to  understanding  physical  mechanisms. 


1.2  Large  Eddy  Simulation 

Large  Eddy  Simulation  is  a  technique  in  which  the  time-dependent  equations  are  solved  for  the 
large  scale  motions.  The  small  scale  turbulence  that  cannot  be  resolved  by  the  mesh  are  approx¬ 
imated  using  model  aptly  referred  to  as  a  “subgrid-scalc  model”.  Such  an  approach  is  grounded 
in  the  realizatino  that  most  of  the  energy  is  carried  by  the  larger  eddies  and  current  computing  ca¬ 
pacity  is  sufficient  to  solve  the  time-dependent  equations  that  can  resolve  these  eddies.  According 
to  Chapman  (22),  the  computational  cost  of  LES  for  wall-bounded  flows  with  near-wall  resolution 
increases  with  the  Reynolds  number  as  Re}  H.  However,  boundary-layer  grids  should  still  suffi¬ 
ciently  resolved  in  the  wall-parallel  directions  to  accurately  resolve  the  near- wall  eddies,  which 
are  in  fact  the  large  eddies  close  to  a  solid  surface.  Notwithstanding  the  use  of  the  best  wall- layer 
treatment,  LES,  is  far  from  affordable  in  complex  aerodynamic  calculations  (49). 

Subgrid  scale  (SOS)  models  are  typically  based  on  the  fact  that  the  global  dissipation  level  is 
set  by  the  largest  eddies  and  the  smaller  eddies  exhibit  similar  behavior  while  transmitting  energy 
to  the  smallest  scales  (25).  Therefore,  the  energy  content  carried  by  smaller  eddies  is  represented 
statistically,  using  a  closure  model.  The  small  scale  quantities  are  separated  from  the  large  scale 
quantities  by  means  of  a  spatial  filtering  operation, 

/(f)  =  I  f(x)G(x  -  x'\  A)dx‘ ,  (9) 


where  G  is  the  filter  function  and  A  is  the  filter  width. 

Applying  the  filtering  operation  to  the  momentum  equations. 


Su,  5  .  _  ,  1  Sp  dr,.  52ut 

- - 1-  —  (W,U,)  —  — “T—  -  7 -  +  -  - 

St  Sxj  3  pSxj  dXj  SxjSxj 


(10) 


where  the  subgrid  scale  stress  ti}  —  u,Uj  —  u,itj.  Traditionally,  the  SGS  stresses  that  are  modeled 
are  divided  into  three  contributions  (25), 
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•  The  resolvable  part  -  Leonard  stresses  L, 3  =  u,Uj  —  u,Uj 


•  The  cross  terms  Cv  —  u,  «' 

•  The  SGS  Reynolds  stresses,  RtJ  =  u'u' 

Hddy  viscosity  models  are  commonly  used  for  closure,  and  these  relate  the  anisotropic  part  of  the 
SGS  stress  tensor  to  the  strain  rate  tensor: 


1.2.1  Subgrid  Scale  Models 

A  number  of  SGS  models  have  been  developed  with  the  most  popular  including  the  Smagorinsky 
model,  scale  similarity  models,  Renormalization  Group  models  (RNG)  theory  models,  and  dy¬ 
namic  sub-grid  scale  models  (20).  One  of  the  commonly  used  subgrid  scale  model  based  on  the 
eddy  viscosity  concept  is  the  Smagorinsky  model.  The  model  lengthscalc  is  given  by  l  =  CsA 
where  Cs  is  the  Smagorinsky  constant  and  A  being  the  filter  width.  The  velocity  scale  is  obtained 
by  assuming  that  the  smaller  scales  are  in  equilibrium  and  the  energy  obtained  from  the  larger 
scales  is  dissipated  instantaneously  and  completely  to  the  smaller  scales. 

=  C'jA2|S|,  \S\  =  (2SijSi})*  (12) 

1.3  Reynolds-averaged  Navier-Stokes 

RANS  techniques  are  based  on  the  Reynolds  decomposition  according  to  which  a  flow  variable  is 
decomposed  into  a  mean  quantity  and  a  fluctuating  quantity, 

f(x,t)  =  f(x)  +  (13) 

where  /(. r.  t)  is  the  flow  variable,  f(x)  represents  the  mean  value  of  the  flow  variable  and  f'(x,  t ) 
represents  the  fluctuation  from  the  mean.  When  this  decomposition  is  applied  to  the  Navier-Stokes 
equations,  the  result  is  an  equation  for  the  mean  quantities  containing  an  extra  term  known  as  the 
Reynolds  stress  tensor  (r  —  —  /m'u').  The  Reynolds  stress  tensor  has  nine  components,  and  the 
term  represents  the  average  flux  of  the  j-momentum  along  the  i-dircction,  which  is  also  equal  to 
the  average  flux  of  the  i-momentum  along  the  j -direction  (21).  This  quantity  represents  the  transfer 
of  momentum  by  the  fluctuating  motion  into  the  mean  flow  and  is  subjected  to  closure  modeling. 
The  flow  properties  in  a  RANS  calculation  are  typically  subjected  to  time  averaging  defined  by, 

1  fl+T 

/(£}=  lim  —  /  /(x'.r)dr.  (14) 

r-oo  1  Jt 
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Closure  modeling  for  the  Reynolds  stress  term  is  usually  accomplished  by  the  use  of  simple  eddy 
viscosity  models.  Most  RANS  models  invoke  the  Boussinesq  approximation  and  therefore  imply 
an  analogy  between  molecular  and  turbulent  transport,  linking  the  turbulent  stress  to  the  mean 
rate  of  strain  using  an  eddy  viscosity.  This  assumes  that  the  turbulent  eddies  in  a  flow  transfer 
momentum  in  much  the  same  fashion  as  molecular  interactions  in  a  gas.  One  of  the  drawbacks  of 
this  analogy  is  that  while  molecular  interactions  occur  at  much  smaller  scales  as  compared  to  the 
length  scales  over  which  flow  properties  are  changing,  turbulent  eddies  have  important  interactions 
at  scales  that  are  comparable  to  the  length  scales  of  variation  of  the  mean  motion  of  the  flow  (24).  In 
spite  of  this  rather  rough  assumption,  the  eddy  viscosity  approach  is  widely,  and  often  successfully, 
applied. 

As  described  by  Rodi  (19),  there  are  three  general  types  of  eddy  viscosity  models: 

•  Zero-equation  models  -  vt  is  specified  directly; 

•  One-equation  models  -  the  length  scale  is  specified  while  the  velocity  scale  is  solved  by  a 
transport  equation; 

•  Two-equation  models  -  transport  equations  are  used  to  solve  for  the  length  scale  and  velocity 
scale. 

Zero-equation  models  are  usually  not  preferred  as  they  fail  to  take  into  account  flow  history  effects 
and  assume  that  turbulence  is  dissipated  where  it  is  generated,  meaning  that  there  is  no  transport 
of  turbulence  (19).  One-equation  and  two-equation  models  are  much  more  commonly  employed 
in  engineering  practice. 

RANS  models  do  not  work  well  in  alt  cases  as  they  are  based  on  empiricism  and  calibrated  in 
a  simple  pool  of  thin  shear  flows.  It  would  be  difficult  to  obtain  accurate  results  for  very  complex 
flows  and  flows  characterised  by  massive  separation  using  RANS  approaches  (49).  A  number  of 
different  RANS  models  have  been  considered  in  the  past  with  no  particular  model  emerging  as  a 
superior  model  for  all  types  of  flows.  As  described  by  Hunt  (68),  in  some  cases  the  duration  of 
a  disturbance  is  smaller  than  the  intrinsic  time  scale  of  the  turbulence  and,  consequently,  an  erro¬ 
neous  turbulence  model  would  have  little  effect  on  the  mean  flow.  However,  the  turbulence  model 
could  result  in  inaccurate  calculations  for  certain  types  of  flows  in  which  the  intrinsic  timescale 
is  smaller  than  the  distortion  timescale,  for  example,  in  some  shear  flows  (49).  There  is  also  a 
general  agreement  that  the  ultimate  potential  of  eddy-viscosity  models  does  not  include  separated 
flows  over  three-dimensional  geometries  (49). 

The  majority  of  modem  models  employed  in  practice  are  based  on  solution  of  transport  equa¬ 
tions  for  either  the  eddy  viscosity  itself  (one-equation  models)  or  equations  for  turbulence  scales 
that  arc  used  to  form  the  eddy  viscosity  (c.g.,  two-equation  models  such  as  K.  —  s).  Some  of  the 
models  which  were  used  for  predictions  in  the  current  work  have  been  summarized  in  the  following 
sections. 
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1.3.1  Spalart-Allmaras  One-Equation  Model 

The  Spalart-Allmaras  (referred  to  as  S-A  throughout  the  document)  model  consists  of  a  transport 
equation  for  the  turbulent  viscosity,  developed  using  empiricism,  dimensional  analysis,  Galilean 
invariance  and  selective  dependence  on  molecular  viscosity  (27).  The  S-A  model  is  local  i.e., 
the  equation  at  one  point  does  not  depend  on  the  solution  at  other  points,  and  is  compatible  with 
grids  of  any  structure  (27).  Turbulence  models  developed  prior  to  S-A,  such  as  the  Cebeci-Smith 
model  (63),  Baldwin-Lomax  model  (64)  and  the  Johnson-King  model  (65)  are  boundary  layer 
models  in  spirit,  physically  treating  the  whole  boundary  layer  as  a  single  tightly-coupled  module, 
which  becomes  incorrect  in  the  presence  of  detached  and  multiple  shear  layers.  In  their  implemen¬ 
tation,  all  of  these  earlier  models  relied  on  velocity  and  vorticity  profiles  varying  along  a  smooth 
gridline,  thus  being  non- local  and  restricted  to  structured  grids.  Further,  these  are  algebraic  mod¬ 
els,  and  there  are  occasions  such  as  the  presence  of  two  solid  bodies  in  the  flowfield  which  would 
require  decisions  that  would  not  facilitate  automation  of  the  method.  Some  models  such  as  the  k-( 
model  (66)  are  local  but  lack  accuracy  in  prediction  of  shock/boundary  layer  interactions  or  separa¬ 
tions  on  smooth  surfaces.  The  k-t  model  may  be  suitable  in  predicting  massively  separated  flows, 
but  require  finer  grids  near  a  wall,  involve  strong  source  terms  that  often  degrade  the  convergence, 
and  demand  non-trivial  upstream  and  freestream  conditions  for  the  turbulence  variables.  To  offset 
the  ncar-wall  problems,  wall  functions  are  typically  introduced,  though  which  lose  any  justification 
in  separated  flows  (27).  While  the  Baldwin-Barth  model  (67)  is  an  attractive  intermediate  which  is 
local,  and  at  the  same  time  does  not  require  a  finer  resolution  than  the  velocity  field  itself. 

The  idea  of  developing  the  S-A  model  was  prompted  by  the  opinion  that  a  one-equation  model 
generated  as  a  simplified  version  of  the  k-e  model  is  not  optimal  and  that  a  one-equation  model 
generated  afresh,  enabling  fuller  control  over  the  performance  of  the  model.  For  example,  the 
Baldwin-Barth  model  was  derived  from  the  k-e  model,  via  assumptions  and  its  diffusion  term  is 
constrained  by  virtue  of  the  assumptions  made. 

The  S-A  model  has  similar  properties  as  that  of  the  Baldwin-Barth  model  in  terms  of  compati¬ 
bility  with  unstructured  grids  and  near-wall  behavior  and  is  more  robust.  There  are  four  versions, 
with  the  simplest  applicable  to  free  shear  flows,  to  the  most  complete  applicable  to  viscous  flows 
past  solid  bodies  and  with  laminar  regions. 

The  S-A  model  solves  an  equation  for  the  variable  0  which  is  dependent  on  the  eddy  viscosity. 
The  model  has  been  derived  empirically  rather  than  using  transport  equations,  and  includes  a  wall 
destruction  term  that  reduces  the  eddy  viscosity  in  the  laminar  sub-layer.  The  model  also  contains 
trip  terms  to  provide  smooth  transition  to  turbulence  if  required.  The  model  is  given  by: 


=  CM  {1  -  fa)  [v-CC^  +  *>)  V  »)  +  Q>2(V^)a]  “ 


r  _  ^bl  r 
CwlJw  2 


■/fiA  U2. 
(15) 


6 


The  eddy  viscosity  is  obtained  via. 


Vt  =  Vjv  1  , 


/til  — 


XJ  +  eti 


where  u  is  the  molecular  viscosity.  The  production  term  is  given  by. 


(16) 


S  =  S 


Kzd 2 


fv2  i  fv 2  1 


1  +  Xfvl 


(17) 


where  5  is  the  magnitude  of  vorticity  and  d  is  the  distance  to  the  closest  wall.  The  wall  destruction 
function  fw  is  given  by. 


1x0=9 


9  =  r  +  cw2(r 6  -  r) , 


v 

Sk2cP  ' 


(18) 


The  trip  functions  in  the  model  are  given  by, 

Jt\  ~  cngtexp  f-c*2 ^5  [d2  +  9?df] ^  ,  gt  =  min  ^0. 1 ,  ( 1 9) 

The  trip  function  /(1  is  specified  in  terms  of  the  distance  dt  from  the  field  point  to  the  trip,  the  wall 
vorticity  at  the  trip  and  A U  which  is  the  difference  between  the  velocity  at  the  field  point  and 
that  at  the  trip.  The  trip  function  ft2  is  given  as, 


ft 2  -  c(3e-c'^2 


(20) 


The  wall  boundary  condition  is  v  =  0.  The  constants  are  cbi  —  0.1355,  a  =  5,  cb2  =  0.622, 
k  =  0.41,  cw  1  =  cbl/K2  +  (1  +  Cb2)/<J,  cw 2  -  0.3,  cw 3  =  2,  cvi  =  7.1,  cn  =  1,  cl2  -  2,  c(3  =  1.1 
and  Cj 4  =  2. 


1.3.2  Spalart-Allmaras  model  with  rotation  correction 


Spalart  and  Shur  (52)  proposed  a  modification  to  the  Spalart-Allmaras  model  to  account  for  stream¬ 
line  curvature  and  system  rotation.  The  approach  involves  second-order  derivatives  of  the  velocity 
field.  The  S-A  model  with  rotation  correction  consists  of  a  modification  to  the  source  term  (Cb\u>i>) 
that  is  multiplied  by  the  rotation  function  fT\,  where  0  is  the  modified  eddy  viscosity.  The  rotation 
function  is  defined  as, 


In  (r\f)  =  (1  +crj) 


2r* 

1  +  r* 


[l  -  ct3  tan  1  (Cr2r )]  -  Cn 


(21) 


where  r*  and  f  arc  nondimensional  quantities  given  by. 


_  S 

(22) 

iiJ 

r  — 

D4 

C  _  n  -  ( 9u,  duj} 

>J  \dxj  dxij 

| 

=  0,5  ( 

du,  ditj\ 
dxj  dxi)  ‘ 

D2  —  ujjj2  4-  Sjj 2 

(23) 
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1.3.3  Shear  Stress  Transport  model 


The  Shear  Stress  Transport  (SST)  model  was  developed  by  Menter  (74)  to  improve  the  accuracy  of 
the  k-Lj  model  in  the  prediction  of  separated  flows.  The  model  combines  k-t  and  k- u  formulations, 
being  identical  to  the  k-u  model  in  the  inner  region  of  the  boundary  layer  (up  to  approximately 
6/2)  and  gradually  changes  to  k-e  model  away  from  the  walls  and  in  the  outer  wake  region.  A 
blending  function  F]  is  introduced  to  bridge  from  k-u  near  the  wall  to  k-t  in  the  freestrcain.  The 
transport  equations  governing  k  and  u>  take  the  form: 


D{nk) 

Dt 


dui 
'  dij 


=  Tij™  -  0pu)k  +  — 


8 

dx} 


8k 

(ti  +  OkVt) 


Dpu)  7 p  _  fhij 

~Dt 


~ - T, 


fit  >J  dxj  dxj 


3uj 


+  2p{\  —  F\)a,„2 


1  8  k  8u 
u  dxj  dxj 


where  is  the  modeled  turbulent  shear  stress.  The  switching  function  Fj  is  given  by. 


(24) 

(25) 


or#]  =  mm  max 


C  Dku  —  tnax 


F]  =  tanh(arg‘l) 

(  (  yfk  500// 


Apa^k 


0,09un/  J  CDkjJy2 


2  '  ^£;10  » 

(jj  ax,  ax, 


(26) 

(27) 

(28) 


The  switching  (blending)  function  F\  is  also  used  to  determine  the  values  of  the  model  constants. 
If  0|  represents  a  generic  constant  in  the  k-u  equations  and  (f> 2  represents  the  same  constant  in 
the  k-t  equations,  then  the  model  constants  employed  in  the  transport  equations  for  k  and  w  arc 
obtained  by, 

0  =  F]<p,+(l  -F,)02  (29) 


The  turbulent  eddy  viscosity  is  determined  by  17  =  k/u>.  The  SST  model  limits  the  turbulent  shear 
stress  to  pa  jit  where  oi  =  0.31.  The  expression  for  the  eddy  viscosity  can  therefore  be  obtained 


from. 


a\k 


t/t  =  — —112 -  (30) 

max(a]W,  i7F2) 

where  fl  is  the  absolute  value  of  vorticity.  The  function  F2  is  included  to  prevent  singular  behavior 
in  the  frecstream  where  Q  goes  to  zero  and  is  given  by. 


F,  =  tanfi(argi2)argi  =  ^7)  <3I> 

The  k-u  model  constants  are  given  by  ok\  =  0.85,  <twi  —  0.5,  0\  =  0.0750,  0‘  =  0.09,  k  —  0.41, 
7,  =  3i/0‘  -  /\f0*,  The  values  of  the  k-t  model  constants  are  07-2  =  1,  2  —  0.856, 

0,  =  0.0828,  .0*  =  0.09,  k  =  0.41,  72  =  A/ff"  - 
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1.4  Detached  Eddy  Simulation 


RANS  approaches  arc  often  acceptable  in  the  thin  shear  layers  where  the  methods  have  been  cali¬ 
brated.  In  other  regimes,  especially  flows  in  which  the  turbulent  eddies  are  not  “standard”,  i.e.,  not 
in  the  calibration  range  of  the  model,  the  performance  of  RANS  models  is,  at  best,  uneven.  This 
in  turn  motivates  other  strategies,  one  of  the  primary  alternatives  being  Large- Eddy  Simulation 
(LES).  The  application  of  LES  to  prediction  of  turbulent  flows  in  practical  configurations  is  in¬ 
creasing,  primarily  in  internal  flows  and  especially  to  flows  that  include  chemistry  or  contain  more 
than  one  phase.  Unfortunately,  the  computational  cost  which  arises  from  the  application  of  LES 
to  a  complete  configuration  such  as  an  airplane,  submarine,  or  road  vehicle  at  practical  Reynolds 
numbers  is  prohibitive.  The  high  cost  of  LES  arises  because  of  the  resolution  requirements  for  the 
boundary  layers,  an  issue  that  remains  even  with  fully  successful  wall-layer  modeling  (49). 

Hybrid  methods  combine  RANS  and  LES  techniques  with  the  aim  of  capitalizing  on  the  rela¬ 
tive  strengths  of  each  approach.  Perhaps  the  most  widespread  RANS-LES  method  in  use  today  is 
Detached- Eddy  Simulation  (DES).  DES  was  originally  proposed  by  Spalart  et  al.  (26)  as  a  cost- 
effective  and  plausibly  accurate  approach  for  predicting  flows  experiencing  massive  separation. 
The  method  was  originally  intended  to  predict  the  entire  boundary  layer  using  a  RANS  model  and 
with  an  LES  treatment  intended  for  the  separated  regions.  DES  has  performed  extraordinarily  well, 
enabling  computationally  feasible  predictions  of  a  range  of  complex  flows  at  high  Reynolds  num¬ 
bers  that  are  either  difficult  to  model  accurately  using  RANS  or  impose  a  computational  cost  that 
prevents  the  use  of  LES  ( 1 6),  (18),  (17),  (57).  Given  the  successful  applications  of  the  technique  to 
date,  there  is  now  strong  interest  in  expanding  the  range  of  applications  which  may  be  accurately 
predicted  using  DES. 

In  the  1997  version  of  the  method,  DES  and  for  that  matter  most  other  hybrid  RANS-LES 
methods,  the  strategy  is  modification  of  the  turbulence  model  to  ensure  a  RANS  treatment  of  the 
(thin)  attached  boundary  layers  and  LES  treatment  of  the  detached  boundary  layers.  The  original 
formulation  of  the  mode!  proposed  by  Spalart  et  al  (26)  achieved  this  behavior  by  modifying  the 
length  scale,  d  (the  wall  distance  in  the  case  of  S-A),  that  enters  the  turbulence  model  and  controls 
the  eddy  viscosity.  The  S-A-based  DES  formulation  is  achieved  by  re -defining  the  length  scale  in 
the  S-A  RANS  model  (shown  here  without  the  trip  terms)  from  d  to  d  in  the  destruction  term  of 
the  model: 


DV 

Dt 


cm  Sis  cu,  i 


V 

2  1 

H - 

jl. 

V((t/  +  i/)V5)  +  C42IVPI' 


If  the  production  and  destruction  terms  of  the  model  are  dominant,  then 


and  consequently, 


C41  Sis  —  cUJ]  ^ 


v  <x  Sd 2 . 


(32) 


(33) 


(34) 
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In  the  Smagorinsky  model  the  eddy  viscosity  is  proportional  to  the  local  rate  of  strain  and  filter 
width,  i.e,,  t/sgs  —  (Cs A)2|S|  and  the  subgrid  length  scale  is  lsgs  —  Cs A  where  Cs  is  the  Smagorin¬ 
sky  constant  and  A  is  the  filter  width,  typically  taken  as  some  measure  of  the  local  grid  spacing. 
Equation  (34)  shows  that  a  Sinagorinsky-like  form  is  achieved  if  d  is  made  proportional  to  the  grid 
spacing.  In  particular,  the  DES  formulation  is  achieved  by  prescribing  the  model  length  scale  as, 

d  =  min(d,  Cdes&)  >  A  =  max(Ax,  Ay,  Az) ,  (35) 

Foes  —  ^  (36) 

where  again  d  is  the  distance  to  the  nearest  wall,  Cdes  is  a  model  constant  to  be  determined,  and 
A  is  a  measure  of  the  local  grid  spacing  taken  as  the  largest  of  the  three  grid  spacings  in  each  of  the 
coordinate  directions.  Note  also  that  the  length  scale  should  be  changed  from  d  to  d  everywhere 
it  appears  in  the  S-A  model.  Near  a  solid  surface  the  wall-parallel  grid  spacings  determine  A  (for 
RANS-type  boundary  layer  grids)  and  are  larger  than  the  distance  to  the  wall  d.  Thus,  within 
the  boundary  where  d  <S  A  boundary  layer  properties  arc  predicted  using  the  S-A  RANS  model. 
Away  from  the  wall  where  Cdes&  <  d,  a  Smagorinsky- like  SGS  model  is  obtained.  The  constant 
Cdes  was  set  by  Shur  et  ai  (56)  in  computations  of  homogeneous  and  isotropic  turbulence.  The 
parameter  Fdes  36  is  defined  as  the  ratio  of  the  turbulent  lengthscale  to  the  wall-distance  and  is 
exactly  equal  to  1  at  the  RANS-LES  interface. 

As  shown  above,  DES  is  based  on  a  single  turbulence  model  and  is  not  zonal,  i.e.,  with  explicit 
declarations  of  zones  as  RANS  or  LES.  The  change  in  the  length  scale  from  the  RANS  form  to 
that  used  in  DES  leads  to  a  model  that  becomes  region-dependent.  In  most  cases  a  RANS  model 
is  active  in  the  entire  boundary  layer  with  a  subgrid-scale  model  active  away  from  solid  surfaces. 
A  key  feature  that  is  common  to  classical  LES  is  incorporation  of  the  grid  spacing,  through  the 
appearance  of  A.  This  is  analogous  to  the  foundation  of  LES  in  which  a  there  is  a  filter  width  that 
controls  the  end  of  the  energy  cascade,  and  can  be  reduced  in  order  to  increase  the  range  of  scales 
and  therefore  improve  the  physics  of  the  simulation  (49). 

The  RANS  and  LES  regions  are  separated  by  an  interface  that  corresponds  to  the  location  where 
the  model  length  scale  becomes  proportional  to  A.  This  interface  is  dictated  by  the  grid  and  can 
result  in  degradation  of  DES  predictions.  In  the  natural  DES  applications  envisioned  by  Spalart  et 
al.  (26),  the  entire  boundary  layer  is  in  “RANS  mode”,  i.e.,  with  a  RANS  model  active,  the  RANS- 
LES  interface  is  outside  the  boundary  layer,  and  there  are  no  adverse  effects  of  the  interface  on 
boundary  layer  predictions.  Grid  refinement  in  the  wall  parallel  directions  (both  streamwise  and 
spanwise)  will  cause  the  RANS-LES  interface  to  move  closer  to  the  wall  and  possibly  within 
the  boundary  layer.  This  can  activate  the  length  scale  switch  to  the  LES  branch  ( d  =  Cdes  A) 
and  reduce  eddy  viscosity  levels  below  that  which  wouid  be  obtained  from  the  RANS  model. 
The  natural  applications  summarized  in  the  next  section  focus  on  flows  that  experience  massive 
separation.  The  flow  around  a  circular  cylinder,  for  example,  is  representative  in  that  a  boundary 
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layer  thinned  by  acceleration  along  the  body  detaches  into  a  region  of  massive  separation.  On 
grids  that  are  typical  for  RANS  prediction  of  the  boundary  layers  in  these  and  other  flows,  errors 
resulting  from  the  RANS-LES  interface  are  not  encountered.  A  major  issue  in  a  DES  is  that  a 
separating  shear  layer  needs  to  generate  random  eddies  which  it  did  not  possess  in  the  boundary 
layer.  This  process  of  generating  three-dimensional  structures  (“LES  content”)  is  more  easily 
accommodated  by  a  thin  shear  layer  rapidly  departing  from  a  wall.  In  massively  separated  flows  at 
high  Reynolds  numbers  instabilities  in  the  detaching  shear  layers  are  rapidly  amplified  and  three- 
dimensional  structures  quickly  fill  the  wake. 

Recently,  Spalart  el  al.  (23)  have  addressed  the  issue  for  natural  DES  applications  where  it 
is  preferable  to  over-ride  the  DES  limiter  and  maintain  RANS  behavior  in  the  boundary  layers, 
independent  of  At|  relative  to  5.  This  ensures  that  the  eddy  viscosity  levels  would  be  given  by  the 
RANS  model  in  flows  with  thick  boundary  layers  and  shallow  separation,  thereby  maintaining  the 
correct  levels  of  the  total  stress  and  skin  friction.  The  issue  addressed  by  Spalart  el  al.  (23)  is 
illustrated  using  Figure  1.  Shown  in  the  figure  are  three  boundary  layer  grids.  Figure  la  shows 
a  “Type  I”  grid  and  is  typical  of  RANS  and  of  DES  with  a  thin  boundary  layer  where  the  wall- 
parallel  spacings  Ax  or  A  z  set  the  length  scale  A,  i.e.,  A  =  maxfAx,  Ay,  Az)  and  the  DES  length 
scale  d  =  d  throughout  the  boundary  layer  since  Cdes&  >  d  f°r  Type  I  grid.  This  in  turn 
implies  that  the  RANS  mode!  is  maintained  throughout  the  boundary  layer  and  DES  functions  as 
intended,  i.e.,  an  LES  treatment  is  bypassed  in  large  areas  of  thin  boundary  layers. 


(a) 

(*) 

i 

(c) 

^ _ ^  : 

_ 

Figure  1 :  Boundary  layer  grids,  (a)  Type  I  grid  for  natural  DES;  ( b )  Type  II  ambiguous  grid;  (c) 
Type  Ill  grid  for  wall-layer  modeling  in  LES. 


The  other  extreme  is  Figure  Ic  which  shows  a  “Type  II”,  grid  with  the  wail-parallel  grid  spac¬ 
ings  much  smaller  than  6  (A  as  5/20  is  a  reasonable  value  for  wall-modeled  LES  (72)).  In  the  Type 
Ml  grid,  the  closure  is  a  subgrid  model  throughout  the  bulk  of  the  boundary  layer  (d  —  C'pesA), 
and  a  RANS-like  wall-layer  model  very  close  to  the  surface.  Between  the  RANS-modeled-rcgion 
and  LES  region  is  a  grey  area  which  has  been  the  subject  of  significant  attention  (72).  Type  III  grids 
are  created  by  design  for  wall-modeled  LES  applications  and  related  studies  have  been  reported  by 
Nikitin  et  al.  (72)  and  Piomelli  et  al.  (40). 


11 


Problems  arise  on  the  “ambiguous"  Type  II  grid  shown  in  Figure  lb.  The  grid  resolution  is 
sufficient  to  activate  the  DES  limiter  (i.e.,  cause  d  =  Cdes&)  within  the  boundary  layer  though 
mesh  spacings  are  not  fine  enough  to  support  resolved  velocity  fluctuations  (i.e.,  LES  content).  In 
this  case  the  DES  limiter  reduces  the  eddy  viscosity,  and  the  modeled  Reynolds  stress,  without  the 
generation  of  any  sizable  resolved  stress  to  restore  the  balance.  Spalart  el  al.  (23)  refer  to  this  as 
“Modeled  Stress  Depletion”,  or  MSD.  This  occurs  when  the  grid  is  gradually  refined  starting  from 
Type  I,  e.g.,  w'hen  the  user  is  seeking  grid  convergence,  or  when  geometry  features  demand  a  fine 
wall-parallel  grid.  The  regime  represented  by  the  Type  II  grid  also  occurs  when  a  boundary  layer 
thickens  and  nears  separation. 

Figure  2  shows  the  effect  of  the  adverse  interactions  of  ambiguous  grids  in  DES  of  the  turbulent 
flow  over  a  flat  plate  with  zero  pressure  gradient.  The  Reynolds  number  of  the  flow  based  on  the 
momentum  thickness  of  the  boundary  layer  at  the  end  of  the  flat  plate  was  approximately  9000 
and  the  flow  conditions  corresponded  to  a  Mach  number  of  0.2.  The  two-dimensional  flow  was 
computed  on  structured  grids  of  different  mesh  densities  using  Cobalt  (36).  The  baseline  grid 
consisted  of  100  x  108  points  in  the  streamwise  and  wall-normal  directions  respectively  and  was 
typical  of  a  Type  I  grid  discussed  above.  The  baseline  grid  was  improved  upon  by  increasing  the 
mesh  density  in  the  streamwise  direction  to  yield  grids  consisting  of  300  x  108,  500  x  108  and 
700  x  108  points.  Additionally,  the  wall-normal  spacing  was  altered  on  one  of  the  grids  to  yield 
the  finest  grid  consisting  of  500  x  200  points.  The  initial  wall-normal  distance  on  all  of  the  grids, 
corresponded  to  an  average  y+  of  less  than  unity. 

Figure  2a  shows  the  location  of  the  RANS-LES  interface  on  each  of  the  grids  in  relation  to  a 
RANS  boundary  layer  and  RANS  eddy  viscosity  distribution  at  the  edge  of  the  plate.  The  parame¬ 
ter  Foes  refers  to  the  ratio  of  minfCoEsA,  d)  to  d  and  assumes  a  value  of  unity  at  the  interface. 
The  interface  is  located  well  within  the  boundary  layer  for  ali  of  the  grids  except  the  baseline  grid 
for  w'liich  it  is  approximately  at  the  edge  of  the  boundary  layer.  This  has  a  marked  effect  on  the 
coefficient  of  skin  friction,  mean  velocity,  and  eddy  viscosity  ratio  which  are  plotted  in  Figure  2b, c 
and  d  respectively.  These  plots  indicate  the  effect  of  the  lengthscale  switch  occurring  well  within 
the  boundary  layer  in  DES,  leading  to  a  drop  in  the  levels  of  the  eddy  viscosity  and  modeled 
stresses  as  compared  to  the  RANS  levels,  thereby  altering  the  growth  of  the  boundary  layer  and 
resulting  in  a  poor  prediction  of  the  skin  friction  coefficient. 

The  issue  of  MSD  exists  in  any  hybrid  RANS-LES  method  such  as  Limited  Numerical  Scales 
(LNS)  (43),  which  also  introduces  the  grid  spacing  into  the  turbulence  model  in  order  to  achieve  an 
LES  treatment.  It  is  to  be  noted  that  two  recently  proposed  methods,  Scale-Adaptive  Simulations 
(SAS)  (42)  and  Turbulence-Resolving  RANS  (TRRANS)  (41)  are  based  purely  on  RANS  mod¬ 
els  and  do  not  incorporate  an  explicit  dependence  of  the  mesh  into  the  turbulence  model  though 
computations  performed  using  these  methods  exhibit  some  LES  properties  (e.g.,  a  range  of  eddies 
captured  on  the  grid).  However,  these  methods  are  not  yet  completely  understood  and  have  not.  in 
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Figure  2:  DES97  of  a  flat  plate  boundary  layer,  (a)  FDbs  along  with  profiles  of  the  RANS  velocity 
profile  and  eddy  viscosity  ratio;  ( b )  skin  friction  coefficient;  (c)  mean  velocity;  (</)  eddy  viscosity 
ratio. 
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the  case  of  SAS,  been  able  to  produce  LES  content  in  channel  flow.  Similarly,  TRRANS  fails  to 
sustain  unsteadiness  over  a  backward  facing  step  (41).  Solutions  to  MSD  have  been  proposed  that 
make  DES  “zonal"  by  disabling  the  DES  limiter  in  selected  regions  where  an  attached  boundary 
layer  is  expected.  This  user  intervention  is  adequate  in  simple  geometries  though  defining  such  re¬ 
gions  in  complex,  three-dimensional  configurations  is  difficult  and  could  be  prone  to  errors.  Other 
solutions  to  MSD  have  attempted  to  use  additional  information  from  the  mesh,  such  as  using  the 
aspect  ratio  of  the  grid  cells  in  order  to  dictate  the  location  where  the  model  switches  from  RANS 
to  LES.  Boundary- layer  grid  cells  are  characterized  by  relatively  large  aspect  ratios  and  in  these 
approaches  an  aspect  ratio  much  larger  than  1  is  used  to  indicate  the  boundary  layer  and  the  length 
scale  d  is  biased  to  maintain  its  RANS  value. 

A  drawback  of  this  approach  is  that  in  the  case  of  unstructured  grids,  there  might  be  abrupt 
transitions  from  high-aspect-ratio  hexahedral  cells  near  the  walls  to  low-aspect-ratio  tetrahedral 
cells  away  from  the  walls.  This  could  in  turn  lead  to  relatively  abrupt  switching  between  RANS 
and  LES  modes  and  possibly  large  changes  in  the  DES  length  scale  across  such  an  interface. 
Conversely,  structured  grids  often  have  regions  with  rather  fine,  high-aspect  ratio  grid  cells  outside 
boundary  layers,  which  in  turn  further  complicates  this  approach. 

Forsythe,  Fremaux  and  Hall  (44)  implemented  a  method  in  which  the  DES  lengthscale,  min{d, 
Cdes A)  was  modified  to  extend  RANS  behavior  for  most  of  the  boundary  layer  when  thicker 
boundary  layers  were  encountered.  This  was  done  by  replacing  the  DES  lengthscale,  min(d. 
Cogs  A)  by  a  function  that  has  the  same  limiting  behaviors  when  ri  <<  Coes  A  and  when 
d  »  Coes  A,  but  follows  rf  when  it  is  only  somewhat  larger  than  Coes  A.  The  specific  form 
of  the  lengthscale  switching  function  was  d  =  min(Coeernax(7i2CoesA;,!/d,  A),f/)  with  n  =  3 
The  shortcoming  of  this  method  is  that  the  use  of  fi/CoesA  to  identify  boundary  layers  is  not 
significantly  robust  as  compared  to  the  use  of  cell  aspect  ratios.  It  would  be  difficult  for  users  to 
anticipate  the  thickness  of  the  boundary  layer  at  the  time  of  grid  generation,  especially  in  flows 
that  are  on  the  verge  of  separation.  Grid-adaption  should  improve  the  situation,  but  the  present 
algorithms  produce  strong  local  variations  of  cell  size  and  aspect  ratio.  Mcnler  and  Kuntz  (37) 
proposed  a  solution  within  the  framework  of  DES  with  the  Fj  or  Fa  functions  of  the  SST  two- 
equation  RANS  model  employed  to  delay  the  switch  to  LES  within  the  boundary  layer.  Spalart  et 
al.  (23)  recently  presented  a  newer  version  of  DES  based  on  the  S-A  RANS  model  which  de¬ 
lays  the  switch  from  RANS  to  LES  within  the  boundary  layer.  This  is  discussed  in  the  following 
section. 
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2  Delayed  Detached  Eddy  Simulation 


The  new  version  of  the  model,  Delayed  Detached  Eddy  Simulation  (DDES),  modifies  the  length 
scale  d  in  order  to  preserve  RANS  treatment  of  the  boundary  layer  irrespective  of  the  grid  spacings. 
Essentially,  the  spirit  of  the  modification  is  to  utilize  information  concerning  the  length  scale  of 
the  turbulence  as  predicted  by  the  model,  in  addition  to  the  wall  distance  and  local  grid  spacing 
to  determine  the  location  of  the  interface.  To  obtain  the  DDES  formulation  the  parameter  r ,{  is 
introduced, 


i/t  4-  v 

Td  =  <ju^urJK2di' 


(37) 


where  £/*, }  are  the  velocity  gradients.  Similar  to  r  in  the  S-A  model,  this  parameter  equals  1  in  a 
logarithmic  layer,  and  falls  to  0  gradually  towards  the  edge  of  the  boundary  layer.  The  addition  of 
v  in  the  numerator  corrects  the  very-near-wall  behavior  by  ensuring  that  rcl  remains  away  from  0. 
The  quantity  is  used  as  an  argument  to  the  function, 


U  =  1  -  tanh([Cr,f]n) , 


(38) 


which  is  designed  to  be  1  in  the  LES  region,  where  rj  <£  1,  and  0  elsewhere  (and  to  be  insensitive 
to  r( i  exceeding  1  very  near  the  wall).  Values  of  8  and  3  were  selected  for  C  and  n  based  on  intuitive 
shape  requirements  for  fd,  and  on  tests  of  DDES  in  the  flat-plate  boundary  layer  (23).  The  DES 
length  scale  d  is  then  re-defined  as, 

d  =  cl  -  fd  max(0,  d  —  CoESA) .  (39) 

Setting  f,t  to  0  yields  RANS  (d  =  d),  white  setting  it  to  1  gives  DES97.  For  DES  based  on  most  of 
the  possible  RANS  models,  DDES  will  consist  in  multiplying  by  /,/  the  term  that  constitutes  the 
difference  between  RANS  and  DES. 

The  two-dimensional  flow  over  a  fiat  plate  with  zero  pressure  gradient  was  computed  using 
DDES  with  the  same  flow  parameters  as  those  used  in  the  DES  study  2.  DES97  and  DDES  pre¬ 
dictions  of  the  flow  over  an  Aerospatiale  A-airfoil  at  an  angle-of-atlack  of  13.3°  were  obtained  on 
different  grids  to  focus  on  the  aspect  of  grid-dependence  of  the  solutions.  These  tests  also  enabled 
to  reinforce  the  suitability  of  the  model  constants  (C  and  n  in  equation  38)  in  DDES.  Based  on  the 
investigations,  the  mean  velocity  profiles  at  locations  normal  to  the  airfoil  surface  are  essentially 
grid-independent  in  the  DDES  predictions,  contrary  to  the  trend  indicated  by  DES97.  Figure  3a 
shows  the  distribution  of  1  -  fd  for  different  values  of  the  coefficients  ((10.  4), (8. 2)  and  (8.  3)) 
along  with  the  destruction  term  in  the  S-A  mode!  for  the  boundary  layer  flow  over  a  flat  plate.  The 
distribution  of  1  —  fd  for  the  case  with  coefficients  8  and  3  indicate  the  suitability  of  the  chosen  val¬ 
ues.  Figure  3b  shows  the  distribution  of  1  -  fd  for  different  values  of  the  coefficients  ((10, 4}, (8, 2) 
and  (8.3))  along  with  the  destruction  term  in  the  S-A  model  for  the  flow  over  the  A-airfoil  at 
Re  =  2  x  lO**  and  at  an  angle-of-attack  ofl3.3°.  Figure  3c  and  d  show  the  variation  in  the  profiles 
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Figure  3:  (a)  Model  parameter  1  —  fd  along  with  the  destruction  term  in  the  S-A  model  for  a  fiat 
plate  boundary  layer;  ( b )  Model  parameter  1  —  /rf  along  with  the  destruction  term  in  the  S-A  model 
for  the  flow  over  an  airfoil  at  Re  =  2  x  106  at  an  angle- of-attack  of  13.3°;  (c)  model  parameter 
1  —  fd  from  the  flow  over  an  airfoil  at  Re  =  2  x  106  at  an  angle  of  attack  of  13.3°  with  n  =  -1  and 
C  varied;  (d)  model  parameter  1  —  fd  from  the  flow  over  an  airfoil  at  Re  =  2  x  106  at  an  angle  of 
attack  of  13.3°  with  C  =  8  and  n  varied, 
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of  1  —  /,(  for  various  values  of  C  and  n  for  the  flow  over  the  A -airfoil  at  Re  =  2x  10{i  and  at  an 
angle-of-attack  of  13.3°.  A  larger  value  of  C  ensures  that  fd  is  equal  to  zero  for  a  longer  extent 
of  the  boundary  layer,  while  a  larger  n  increases  the  sharpness  of  the  function,  ensuring  that  fl{ 
quickly  changes  from  1  to  zero  at  the  interface.  Figure  4a  and  b  show  contours  of  the  model  para¬ 
meter  1  —  fd  for  the  flow  over  a  cylinder  at  Re  =  8  x  106  with  C  =  8  and  n  =  2  and  3.  It  can  be 
seen  from  the  contour  plot  that  with  C  =  8  and  n  =  2,  the  transition  region  where  fd  changes  from 
0  to  1  is  wider  as  compared  to  the  case  with  C  =  8  and  n  =  3.  Figure  4c  and  d  show  the  same 
effect  when  n  is  varied  for  the  flow  over  the  A-airfoil  at  Re  =  2  x  106  and  at  an  angle-of-attack  of 
13.3°. 


Figure  4:  Contours  of  the  model  parameter  1  -  /,/  from  the  flow  over  a  cylinder  at  ffe-8x  106: 
(a)  C  =  8,  n  =  2;  (i>)  C  =  8,  n  =  3;  Contours  of  the  model  parameter  1  -  fd  from  the  flow  over 
an  airfoil  at  Re  =  2  x  106  at  an  angle  of  attack  of  13°:  (c)  C  =  8,  n  =  2;  (d)  C  =  8,  n  =  3. 


Figure  5  shows  distributions  of  the  skin  friction  coefficient  and  the  mean  velocity  at  the  end  of 
the  plate.  The  results  show  grid  independence  and  compare  well  with  the  RANS  on  all  the  grids 
except  for  the  700  x  108  grid  for  which  there  is  only  a  slight  mismatch  with  the  RANS.  Figure  6 
shows  profiles  of  the  parameters  Foes  and  1  “  fd  in  relation  to  a  RANS  boundary  layer  streamwisc 
velocity  and  RANS  eddy  viscosity  distribution  at  the  end  of  the  flat  plate.  The  parameter  F oes  is 
the  ratio  of  the  turbulent  length  scale  in  DDES  to  the  wall  distance^/). 

3  Flow  over  a  Circular  Cylinder  using  DDES 

An  application  of  DDES  to  the  flow  over  a  circular  cylinder  provides  an  opportunity  to  assess 
the  new  method  in  a  flow  for  which  DES97  is  a  very  powerful  and  accurate  approach.  While 
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Figure  5:  DDES  of  a  flat  plate  boundary  layer,  (a)  mean  velocity;  ( b )  skin  friction  coefficient. 


Figure  6;  Turbulent  model  properties  along  with  RANS  mean  streamwise  velocity  and  eddy  vis¬ 
cosity  ratio,  (a)  Foes',  (t)  1  -  Id- 
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the  cylinder  is  not  a  candidate  that  should  be  prone  to  MSD  (boundary  layers  prior  to  separation 
are  thin  and  mesh  spacings  on  typical  grids  will  not  locate  the  RANS-LES  interface  within  the 
boundary  layer),  predictions  should  not  be  degraded  using  DDES  compared  to  those  using  DES97. 

Computations  are  presented  at  a  Reynolds  number  based  on  the  freestream  velocity  and  cylin¬ 
der  diameter  D  of  1.4  x  10s  and  8  x  106.  The  lower  Re  predictions  are  assessed  against  previous 
simulations  while  experimental  measurements  were  used  to  evaluate  the  higher  Re  predictions,  For 
Re  ~  1 .4  x  ID5,  the  DES  predictions  obtained  by  Travin  et  al.  (57)  and  Hansen  and  Forsythe  (28) 
were  available  for  assessment  and  comparison.  Mesh  refinement  was  investigated  at  a  higher 
Reynolds  number.  Re  =  8  x  10s,  the  results  of  which  were  compared  with  experimental  measure¬ 
ments  (70),  (71).  The  compressible  Navier-Stokes  equations  were  solved  on  unstructured  grids 
using  Cobalt.  The  spanwisc  coordinate  of  the  domain  along  which  periodic  boundary  conditions 
were  applied  was  AD  for  both  Reynolds  numbers.  The  calculations  at  Re  =  1.4  x  10s  were  per¬ 
formed  using  a  mesh  that  consisted  of  1.434  x  106  cells  (the  same  grid  as  that  used  by  Hansen  and 
Forsythe  (28))  and  characterized  by  a  spanwisc  grid  spacing  of  A  JD  =  0.10.  The  calculations  at 
Re  =  1.4  x  105  were  used  to  provide  an  initial  assessment  of  DDES  against  DES97  results.  The 
influence  of  mesh  refinement  was  investigated  at  the  higher  Reynolds  number  using  three  grids 
having  1.46  x  106,  3.71  x  106,  and  9.83  x  106  cells  and  are  referred  to  as  “coarse  grid”  or  “eg”, 
“medium  grid"  or  “mg”,  and  “fine  grid”  or  “fg”.  The  medium  grid  was  created  by  refining  in  all 
directions  the  coarse  mesh  by  a  factor  of  \/2.  The  finer  grid  had  another  factor  of  \[2  refinement 
in  each  coordinate  direction.  The  corresponding  spanwise  grid  spacings,  which  determined  A  in 
DES97,  were  A ,/D  =  0.10,  A Z/D  =  0.071,  and  A Z/D  —  0.05  for  the  coarse,  medium,  and 
fine  grids,  respectively.  Figures  7a  and  b  show  the  cross-sections  of  the  mesh  in  the  vicinity  of  the 
cylinder.  The  grid  was  created  using  VGRIDns  (55)  and  is  comprised  of  prisms  near  the  cylinder 
surface  and  tetrahedra  away  from  the  wall.  For  all  of  the  grids,  the  initial  wall  normal  distance 
was  within  one  viscous  unit  on  average  and  a  stretching  rate  of  1.2  was  applied  to  the  wall-normal 
grid  within  the  boundary  layer.  The  domain  extended  approximately  20  cylinder  diameters  from 
the  cylinder  surface.  For  both  Reynolds  numbers,  the  super-critical  flow  is  approximated  by  com¬ 
puting  fully-turbulcnt  solutions.  At  the  inlet  of  the  computational  domain  a  small  level  of  eddy 
viscosity  (corresponding  to  \  =  3)  is  prescribed  that  is  sufficient  to  ignite  the  turbulence  model 
as  the  fluid  enters  the  boundary  layers.  For  all  of  the  simulations  reported  here,  the  time  step  was 
fixed  al  O.OID/Uoo  where  f/w  is  the  freestream  velocity. 

3.1  Predictions  at  Re  =  1.4  x  105 

Shown  in  table  1  are  averaged  quantities  from  DES97  and  DDES  predictions  of  the  flow  at  Re  — 
1.4  x  105:  the  drag  coefficient  Cj,  rms  lift  coefficient  C|.rms,  shedding  Strouhal  number  St,  base 
pressure  coefficient  Cp  and  separation  angle  0scp.  Also  included  are  the  DES97  predictions  of 
Travin  el  al.  (57)  and  Hansen  and  Forsythe  (28)  along  with  the  experimental  measurements 
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reported  in  Roshko  (70).  The  spanwise  grid  spacing  for  the  current  computations  and  those  of 
Hansen  and  Forsythe  (28)  is  10  nodes  per  cylinder  diameter  while  the  structured  mesh  for  the 
fully-turbulent  solution  summarized  in  the  table  from  Travin  et  al.  (57)  employed  a  spanwise 
resolution  of  15  grid  points  per  diameter.  In  general,  the  table  shows  that  there  is  good  agreement 
between  the  present  DES97  and  DDES  predictions  and  the  corresponding  results  from  Travin  et 
al.  (57)  and  Hansen  and  Forsythe  (28).  The  drag  coefficient  from  all  of  the  simulations  is 
on  the  low  end  of  the  measurements  summarized  by  Roshko  (70),  that  correspond  to  a  higher 
Reynolds  number  range  (3.5  x  106  <  Re  <  8.4  x  106).  Table  1  also  shows  that  the  back  pressure 
is  slightly  more  negative  for  the  DDES  which  is  not  inconsistent  with  the  slightly  higher  rms  lift 
that  is  predicted  compared  to  the  DES97  results  (and  also  Travin  et  al.  (57),  which  have  the  same 
rms  lift  as  in  the  current  DES97),  An  additional  contributor  to  the  differences  in  table  1  is  the 
difference  in  sampling  periods. 


Figure  7:  Cross-section  of  the  grids  in  the  vicinity  of  the  cylinder,  (a)  coarse  (upper  half  plane) 
and  medium  (lower  half  plane)  grids;  ( b )  medium  (upper  half  plane)  and  fine  (lower  half  plane). 


Figure  8:  Time  history  of  the  drag  (in  a)  and  lift  (in  b)  force  coefficients,  Re  =  1.4  x  10'. 


Time  histories  of  the  drag  and  lift  coefficients  for  each  model  are  plotted  in  Figure  8a  and  b, 
respectively.  Each  figure  shows  comparable  behavior  for  the  two  models,  including  the  significant 
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Figure  9:  Coefficient  of  pressure  and  skin  friction  coefficient  on  the  cylinder  surface,  Re  =  1 .4  x 

10s. 


cd 

CltTTns 

St 

&scp 

DES97 

0.58 

0.08 

0.29 

0.64 

98° 

DDES 

0.60 

0.11 

0.28 

0.69 

99° 

Travin  el  al. 

0.57 

0.08 

0.30 

0.65 

99° 

Hansen  et  al. 

0.59 

- 

0.29 

0.72 

- 

Roshko 

0.62-0.74 

- 

0.27 

- 

- 

Table  1:  Lift  and  Drag  co-efficients  for  the  cylinder  computations  at  Re  —  1.4  x  10s 
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modulation  in  the  lift  force  as  also  reported  in  Travin  el  al.  (57)  and  Hansen  and  Forsythe  (28). 
These  modulations  require  very  long  sample  periods  for  convergent  statistics,  probably  longer 
than  those  shown  in  figure  8,  which  again  contributes  to  the  differences  in  statistics  summarized  in 
table  I. 


(c)  (d) 

Figure  10:  Contours  of  the  instantaneous  vorticity  from  the  fine-grid  solutions,  Re  —  8  x  106,  (a) 
DES97,  (b)  DDES;  Contours  of  the  instantaneous  eddy  viscosity  ratio  from  the  fine-grid  solutions. 
Re  =  8  x  106,  (c)  DES97  (cl)  DDES. 

Figure  9  shows  good  agreement  of  the  pressure  coefficient  from  the  two  models  along  with  the 
results  of  Travin  et  al.  (57)  and  experimental  measurements  form  Roshko  (70)  and  van  Nunen  (71), 
The  experiments  were  performed  at  a  higher  Reynolds  number  (8.4  x  1 06  in  Roshko  (70)  and 
7.6  x  10u  in  van  Nunen  (71)  than  the  computations.  However,  the  fully  turbulent  treatment  of 
the  flow  can  be  considered  a  model  of  the  high  Reynolds  number  experiment  so  long  as  transition 
occurs  far  upstream  of  separation.  The  Cf  plot  indicates  that  the  simulation  results  are  far  from 
the  measurements  prior  to  separation.  While  the  simulations  are  of  a  fully  turbulent  boundary 
layer  prior  to  separation,  the  experiments  indicate  the  boundary  layer  was  laminar  even  close  to 
the  separation  line, 
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Figure  1 1:  Coefficient  of  pressure  on  the  cylinder  surface,  Re  =  8  x  106,  (a)  DES97,  (fr)  DDES 

3.2  Predictions  at  Re  =  8  x  106 

Contours  of  the  instantaneous  vorticity  magnitude  in  three  planes  along  the  span  of  the  cylinder 
are  shown  in  Figure  10.  The  visualizations  shown  in  the  figure  are  from  the  fine  grid  DHS97 
and  DDES  predictions.  Characteristic  of  a  natural  DES  application  in  a  flow  experiencing  massive 
separation,  the  figure  shows  that  the  eddy  content  in  the  wake  develops  rapidly  following  boundary 
layer  detachment.  Contours  of  the  instantaneous  eddy  viscosity  ratio,  u{jv  in  a  single  spanwise 
plane  of  the  cylinder  arc  also  shown  in  Figure  10.  Both  DES97  and  DDES  show  similar  levels  of 
eddy  viscosity  in  the  wake. 

The  pressure  coefficient  for  both  models  and  each  grid  is  compared  to  the  experimental  mea¬ 
surements  of  Roshko  and  van  Nunen  in  Figure  1 1 .  Comparison  of  the  frames  shows  that  the  models 
predict  similar  pressure  distributions,  consistent  with  the  close  agreement  observed  in  the  drag  co¬ 
efficient.  The  figure  shows  the  pressure  recovery  in  the  aft  region  for  both  DES  models  is  in  closer 
agreement  to  the  measurements  of  van  Nunen  and  that  the  computations  exhibit  a  lower  minimum 
in  the  vicinity  of  90  degrees  than  measured  in  either  experiment. 

Time  histories  of  the  drag  and  lift  force  coefficients  on  all  three  grids,  averaged  over  the  span- 
wise  dimension,  are  shown  in  Figure  12.  The  histories  shown  in  the  figure  were  subsequently 
averaged  over  80  time  units  to  obtain  the  values  summarized  in  table  2.  The  lift  coefficient  exhibits 
very  large  modulations  on  the  coarse  grid,  which  are  not  observed  on  the  finer  meshes.  The  large 
modulations  on  the  coarse  grid  are  indicative  of  very  weak  thrcc-diincnsionality  due  to  inadequate 
resolution.  The  lift  flatness  factor,  which  provides  a  measure  of  the  irregularity  of  the  signal  with 
respect  to  the  mean  (larger  excursions  from  the  mean  imply  a  larger  lift  flatness  factor)  is  2.21, 
2.16  and  2.04  for  the  coarse,  medium,  and  fine  grids,  respectively.  These  values  are  in  reasonable 
agreement  with  the  value  of  2.1  reported  by  Travin  et  al.  (57)  in  DES  of  the  flow  around  a  cylinder 
at  Re  =  3  x  106. 
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Figure  12:  Time  histories  of  the  drag  and  lift  coefficients  for  DDES  of  the  cylinder,  Re  =  8x  IQ6. 


The  averaged  drag  coefficient,  shedding  Strouhal  number,  and  separation  angle  for  each  of  the 
grids  and  both  DES97  and  DDES  models  are  summarized  in  Table  2.  With  the  exception  of  the 
DDES  prediction  using  the  coarsest  grid,  the  averaged  properties  predicted  by  the  two  models  arc 
in  quite  good  agreement.  In  general,  the  drag  coefficient  is  around  0.4,  which  agrees  well  with  the 
value  of  0,41  reported  in  Travin  el  al.  (57)  from  DES97  predictions  at  Re  =  3  x  10f).  And  with 
the  exception  of  the  coarsest  grid,  the  shedding  Strouhal  numbers  are  0.36  —  0.37,  also  close  to  the 
value  of  0.35  in  Travin  et  al.  (57),  again  at  Re  —  3  x  10^.  The  table  also  shows  that  on  the  finer 
meshes,  boundary  layer  separation  occurs  at  114°,  slightly  more  aft  than  the  value  of  1 11°  in  the 
tolly  turbulent  DES97  prediction  at  Re  —  3  x  106. 


DES97 

St 

Coarse  grid 

0.43 

0.34 

114° 

Medium  grid 

0.40 

0.36 

114° 

Fine  grid 

0.37 

0.37 

1 14° 

DDES 

St 

$  sep 

Coarse  grid 

0.46 

0.33 

115° 

Medium  grid 

0,40 

0.37 

114° 

Fine  grid 

0.38 

0.37 

114° 

Table  2:  Drag  co-efficient,  Strouhal  number  and  angle  of  separation  for  the  cylinder  computations 
at  Re  =  8  x  106 


Figure  13  shows  the  skin  friction  coefficient  over  the  surface  of  the  cylinder  for  the  DDES 
cases.  The  skin  friction  predictions  in  all  cases  show  good  agreement  with  each  other  but  are  far 
above  the  simulation  results  of  Travin  et  al.  (57)  at  Re  =  3  x  106.  Figure  13  also  shows  the 
streamwise  velocity  distribution  in  the  wake  of  the  cylinder  with  a  crinkle-cut  view  of  the  grid.  For 


24 


a  dimensionless  time  step  of  0.01,  the  figure  shows  CFL  numbers  of  0(1)  at  different  locations  in 
the  wake  of  the  cylinder. 


Figure  13:  Time  histories  of  the  drag  and  lift  coefficients  for  DDES  of  the  cylinder.  Re  =  8  x  106. 


Properties  of  the  turbulence  model  are  shown  in  Figure  14.  Plotted  in  figure  14(a)  and  (b) 
are  radial  profiles  of  the  mean  velocity,  eddy  viscosity  ratio,  and  d/d  at  0  =  75°  and  0  =  90°. 
The  quantities  shown  in  the  figure  are  DDES  predictions  for  each  of  the  grids.  Note  that  the 
mean  velocities  are  averaged,  the  profiles  of  the  eddy  viscosity  and  d/d  are  not  averaged.  Though 
not  shown  here,  the  results  from  the  DES97  cases  are  virtually  identical  For  both  locations, 
Figures  14(a)  and  (6)  show  that  the  RANS-LES  interface  (the  location  where  d  becomes  less  than 
d)  moves  closer  to  the  wall  as  the  grid  is  refined,  though  for  this  natural  application  -  thin  boundary 
layers  into  a  region  of  massive  separation  -  the  interface  remains  well  outside  of  the  boundary 
layer.  Figures  14(c)  and  (d)  show  the  radial  profiles  of  the  mean  velocity,  eddy  viscosity  ratio, 
and  1  —  fd  which  is  designed  to  go  to  zero  in  a  region  where  LES  behavior  is  desired.  When 
massive  separation  occurs,  fd  quickly  rises  from  zero  very  close  to  the  wall  enabling  the  growth 
of  instabilities  in  the  LES  region.  Figure  14(e)  and  (/)  show  the  distribution  of  rd  along  with  the 
radial  mean  velocity  profiles  and  eddy  viscosity  profiles.  The  behavior  of  rd  is  based  upon  the 
behavior  of  fd,  quickly  dropping  to  zero  well- with  in  the  boundary  layer  for  all  three  grids. 
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8=90’  8=75" 


Figure  14:  Turbulence  model  properties:  d/d,  (n)  9  =  75°,  (6)  6  =  90°;  Turbulence  model  prop¬ 
erties:  1  -  /,/,  (c)  9  =  75°,  (d)  0  =  90°;  Turbulence  model  properties:  r*,  (e)  9  =  75°,  (/) 
9  =  90° 
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4  Flow  over  the  Aerospatiale  A-airfoil 


The  predictions  of  the  flow  around  the  Aerospatiale  A-airfoil  at  a  chord-based  Reynolds  number 
of  2  x  106  are  presented  in  this  manuscript.  The  flow  at  an  angle-of-altaek  of  a  =  13.3°  is  the 
major  focus  of  the  investigation,  with  a  few  predictions  of  the  flow  at  angles-of-attack  cv  —  7.2° 
and  q  =  12.2°  also  included. 

4.1  Predictions  at  a  =  13.3° 

The  RANS  predictions  of  the  two-dimensional  flow  were  computed  using  the  Spalart-Allmaras 
(S-A)  (with  and  without  rotation  correction)  and  Shear  Stress  Transport  (SST)  turbulence  models. 
Solutions  of  the  fully-turbulent  flow  were  computed  using  both  RANS  models  in  addition  to  pre¬ 
diction  of  the  flow  with  laminar  separation  near  the  leading  edge  using  the  triplcss  approach  within 
the  Spalart-Allmaras  model.  For  the  tripless  case,  transition  from  laminar-to-turbulcnt  flow  along 
the  lower  (pressure)  surface  was  fixed  by  tripping  the  boundary  layer  at  x/C  =■  0.3.  The  pressure 
coefficient  in  the  RANS  solutions  are  in  reasonable  agreement  with  measurements  prior  to  trail¬ 
ing  edge  separation.  The  tripless  solution  cases  using  the  Spalart-Allmaras  model  predict  laminar 
separation  close  to  x/C  —  0.10  with  turbulent  reattachment  close  to  x/C  —  0.12  and  followed 
by  separation  of  the  turbulent  boundary  layer  near  the  trailing  edge  close  to  x/C  ~  0.90.  Trail¬ 
ing  edge  separation  in  the  fully  turbulent  solutions  is  invariably  closer  to  the  experimental  value  of 
x/C  =  0.83.  Comparison  of  the  mean  velocity  predictions  to  measurements  shows  a  behavior  con¬ 
sistent  with  the  mismatch  in  separation  location  compared  to  measurements  in  the  tripless  cases. 
The  rotation-correction  does  not  improve  the  results  appreciably  and  produces  similar  results  as 
the  case  without  the  correction.  DES97  and  DDES  predictions  were  obtained  on  different  grids 
to  focus  on  the  aspect  of  grid-depcndcncc  of  the  solutions.  Based  on  the  investigations,  the  mean 
velocity  profiles  at  locations  normal  to  the  airfoil  surface  are  essentially  grid-independent  in  the 
DDES  predictions,  contrary  to  the  trend  indicated  by  DES97. 

4.2  Approach 

The  flow  over  the  Aerospatiale  A-airfoil  at  an  angle-of-attack  of  13.3°  and  at  a  chord-based 
Reynolds  number  of  2.1  x  106  was  studied  in  great  detail  through  the  LESFOIL  project  (31). 
The  main  advantage  of  using  LES  for  such  a  case  include  the  reduced  sensitivity  to  modeling 
errors  occuring  in  RANS  methods  and  the  use  of  grid  refinement  to  obtain  richer  flow  physics. 
One  of  the  major  findings  of  the  LESFOIL  project  (31)  is  that  even  a  narrow  section  of  the  air¬ 
foil  (at  a  moderately  high  Reynolds  number  and  unswept),  poses  an  extreme  challenge  to  LES  in 
terms  of  computational  cost.  It  was  also  noted  that  other  aspects  such  as  the  numerical  method, 
grid  generation  and  mesh  resolution  requirements,  especially  in  the  near-wall  region,  coupled  with 
considerations  of  the  spanwise  period  in  three-dimensional  eddy-resolving  computations  stress  the 
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capabilities  of  any  simulation  strategy.  The  general  characteristics  of  this  flow  are  illustrated  in 
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Figure  15:  Flow  regimes  on  the  Aerospatiale  A-airfoi!  (31) 


figure  1 5.  As  the  flow  moves  along  the  upper  (suction)  surface  of  the  airfoil,  the  laminar  boundary 
layer  experiences  a  free  transition  in  a  small  laminar  separation  bubble  with  turbulent  reattach  men  I 
at  around  x/C  =  0.12.  The  turbulent  boundary  layer  grows  along  the  suction  surface  of  the  airfoil, 
separating  at  about  x/C  =  0.83,  resulting  in  a  shallow  separation  bubble  that  extends  to  the  trailing 
edge.  Measurements  indicate  a  reverse-flow  region  extending  to  roughly  a  wall-normal  location 
yfC  =  0.016  where  y  defines  the  coordinate  normal  to  the  suction  surface  and  that  the  location  of 
maximum  reverse  flow  occurs  at  y/C  =  0.003.  Along  the  lower  (pressure)  surface  of  the  airfoil, 
the  boundary  layer  was  tripped  at  x/C  —  0.3,  enabling  transition  from  lam inar-to-turbu lent  flow. 
The  measurements  acquired  included  the  airfoil  lift  and  drag,  skin  friction,  pressure  distributions 
as  well  as  the  mean  velocity  and  Reynolds  stress  profiles  at  several  stations  along  the  airfoil  and  in 
the  near-wake  (29),  (30). 

4.2.1  Flow  solver  and  grids 

The  compressible  Navier-Stokes  equations  arc  solved  on  unstructured  grids  using  Cobalt.  The 
numerical  method  used  consists  of  a  cell-centered  finite  volume  approach  applicable  to  arbitrary 
cell  topologies  (36).  The  spatial  operator  uses  the  exact  Riemann  solver  of  Gottleib  and  Groth  (60), 
least-squares  gradient  calculations  using  QR  factorization  to  provide  second  order  accuracy  in 
space,  and  TVD  flux  limiters  to  limit  extremes  at  cell  faces.  A  point  implicit  method  using  analytic 
first-order  inviscid  and  viscous  Jacobians  is  used  for  advancement  of  the  discretized  system.  For 
time-accurate  computations,  a  Newton  sub-iteration  scheme  is  employed,  the  method  is  second- 
order  accurate  in  time.  The  domain  decomposition  library  ParMETIS  (73)  is  used  for  parallel 
implementation  and  provides  optimal  load  balancing  with  a  minimal  surface  interface  between 
zones.  Communication  between  processors  is  achieved  using  Message  Passing  Interface.  The 
working  of  the  solver  has  been  explained  in  greater  detail  in  the  Appendix  section. 
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The  grids  used  in  the  study  were  structured  and  C-type  grids,  generated  using  Gridgen  (32), 
The  grids  are  referenced  as  'coarse  grid  1 ‘coarse  grid  2’,  ‘medium  grid’  and  'fine  grid’  in  table  3. 
For  ‘coarse  grid  1  the  ratio  of  cell  size  in  the  streamwise  direction  to  the  RANS  boundary  layer 
thickness  at  xjC  =  0.5  is  approximately  0.68,  for  the  ‘fine  grid’,  the  corresponding  value  is 
approximately  0.61.  These  values  serve  as  average  measures  of  the  grid  density,  indicating  the 
absence  of  a  marked  difference  between  ‘coarse  grid  1*  and  ‘fine  grid’.  It  must  be  noted  that 
‘coarse  grid  2’  was  formed  by  extruding  ‘coarse  grid  1  ’  in  the  spanwise  direction  with  A,  =  0.005 
and  a  spanwise  extent  of  0.01,  a  fictitious  value  that  does  not  have  a  marked  effect  on  altering 
the  location  of  the  RANS-LES  interface  in  DES97.  Figure  1 6  shows  ‘coarse  grid  1 1  close  to  the 
surface  of  the  airfoil  and  the  spanwise  extent  of ‘coarse  grid  2’. 


Figure  16:  (a)  Coarse  grid  1 ;  (b)  Spanwise  extent  of  coarse  grid  2 
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grid 

2D/3D 

A* 

number  of  cells 

coarse  grid  1 

2D 

- 

76,000 

coarse  grid  2 

3D 

0.005 

76, 000  x  2 

medium  grid 

2D 

85,000 

fine  grid 

2D 

- 

98,000 

Table  3:  Grid  sizes  used  to  study  the  flow  over  the  Aerospatiale  A-airfoil  at  an  angle-of-attack  of 
13.3°  and  at  a  chord-based  Reynolds  number  of  2.1  x  106 


4.2.2  Results 

Measurements  show  that  the  flow  over  the  Aerospatiale  A-airfoil  experiences  a  laminar  separation 
in  the  vicinity  of  the  leading  edge  region,  just  downstream  of  the  peak  negative  pressure  along  the 
suction  side.  Transition  occurs  in  the  separated  shear  layer  with  the  reattached  turbulent  boundary 
layer  evolving  further  along  the  suction  side  prior  to  a  subsequent  separation  near  the  trailing 
edge.  The  laminar  separation  and  transition  is  accounted  for  using  the  tripless  approach  outlined 
by  Travin  (57).  The  tripless  approach  provides  a  means  to  accommodate  the  laminar  separation 
and  transition  in  the  separated  shear  layer,  in  the  present  calculations  represented  by  an  activation 
of  the  turbulence  model.  The  eddy  viscosity  upstream  of  the  airfoil  is  zero,  non-zero  values  are 
seeded  into  the  suction  side  of  the  airfoil  using  a  boundary  layer  trip. 

4.2.3  Tripless  and  Fully  Turbulent  solutions 

The  suction  side  trip  was  location  was  varied  to  be  at  i /C  =  0.08,  x/C  =  0.12,  x/C  —  0.13  or 
x/C  —  0.20.  A  second  boundary  layer  trip  is  located  at  xjC  =  0.3  along  the  pressure  surface,  the 
same  location  as  used  in  the  experiments.  Predictions  of  the  lully-turbulent  flow  are  also  obtained 
using  the  S-A  and  SST  turbulence  models.  The  fully  turbulent  solutions  are  computed  by  seeding 
a  small  level  of  eddy  viscosity  into  the  upstream  flow,  sufficient  to  ignite  the  turbulence  model  as 
the  fluid  enters  the  boundary  layers.  The  initial  (laminar)  separation  near  the  leading  edge  does 
not  occur  in  the  fully  turbulent  solutions.  Other  aspects  such  as  the  results  close  to  the  trailing 
edge  of  the  airfoil  are  also  different  in  the  fully-turbulent  solutions  as  compared  to  the  tripless 
solutions.  In  the  tripless  runs,  there  is  laminar  separation  followed  by  turbulent  rcaltachmenl  close 
to  x/C  =  0.10.  A  second  separation  is  predicted  near  the  trailing  edge  of  the  airfoil.  The  RANS 
separation  locations  are  shown  in  table  4.  The  separation  locations  in  the  fully  turbulent  runs  are 
closer  to  the  experimental  value  o f  x/C  -  0.83.  Among  the  fully  turbulent  runs,  the  SST  model 
predicts  a  value  of  0.82,  closest  to  the  experimental  results,  while  S-A  and  S-A  RC  models  predict 
earlier  separation  (at  x/C  =  0.79  and  x/C  =  0.77  respectively). 

While  the  fully  turbulent  runs  predict  the  location  of  turbulent  separation  to  be  dose  to  the 
experimental  value  of  x/C  =  0.83,  they  fail  to  shed  light  on  the  laminar  separation  that  occurs 
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grid 

RANS 

model 

Top  Trip 

(x/C) 

Transition 

(x/C) 

Turbulent 
Separation  (x/C) 

medium  grid 

S-A  FT 

- 

- 

0.79 

S-A  RC  FT 

- 

- 

0.77 

SST  FT 

- 

- 

0.82 

S-A  tripless 

0.2 

0.106-0.1 18 

0.91 

S-A  RC  tripless 

0.2 

0.102-0.118 

0.89 

fine  grid 

S-A  FT 

- 

- 

0.79 

S-A  RC  FT 

* 

- 

0.77 

SST  FT 

* 

- 

0.82 

S-A  tripless 

0.2 

0.093-0.128 

0,9! 

S-A  RC  tripless 

0,2 

0.102-0.1 18 

0.90 

S-A  tripless 

0.08 

- 

0.88 

S-A  tripless 

0.13 

0.10-0.125 

0.91 

S-A  tripless 

0  C 

- 

0.91 

coarse  grid  1 

S-A  tripless 

0.08 

- 

0.89 

Table  4:  Summary  of  the  RANS  computations  for  a  —  13.3° 


on  the  forward  part  of  the  suction  surface,  motivating  the  need  to  employ  trips  on  the  suction 
surface.  Preliminary  investigations  involved  the  placement  of  the  top  trip  at  a  location  correspond¬ 
ing  to  x/C  ~  0.2.  In  these  cases,  during  the  initial  stages  of  the  simulations,  the  eddy  viscosity 
drifted  upstream  from  the  location  of  the  trip,  inducing  transition  at  a  location  approximately  cor¬ 
responding  to  x/C  ~  0. 1 1.  More  runs  were  conducted  with  the  location  of  the  top  trip  changed 
to  x/C  —  0.13,  x/C  —  0.12  and  x/C  =  0.08.  Application  of  the  top  trip  at  x/C  =  0.13  gave 
a  similar  transition  location,  with  no  laminar  separation  bubble  present  when  the  top  trip  was  lo¬ 
cated  at  x/C  -  0.12  and  x/C  —  0.08.  The  main  shortcoming  of  the  tripless  runs  is  that  they 
predict  a  location  of  turbulent  separation  further  downstream  of  the  experiment  and  of  the  fully 
turbulent  runs.  The  separation  locations  obtained  using  DES97  and  DDES  are  compared  with  S-A 
RANS  predictions  in  table  5,  with  the  location  of  the  top  trip  fixed  at  x/C  =  0.08.  The  DES97 
predictions  show  a  large  variation  in  the  results  using  the  fine  grid  and  either  of  the  coarse  grids. 
The  results  obtained  with  ‘coarse  grid  1  ’  and  ‘coarse  grid  2’  hardly  show  any  variation  indicating 
that  the  spanwise  extent  of ‘coarse  grid  V  hardly  changes  the  location  of  the  RANS-LES  interface, 
the  location  of  which  alters  the  separation  point.  The  DDES  predictions  on  all  grids  arc  in  good 
agreement  among  themselves  and  very  close  to  the  S-A  RANS  results. 

The  pressure  coefficient  over  the  airfoil  on  both  the  pressure  and  suction  sides  using  RANS  is 
shown  in  figure  1 7.  The  figures  show  that  compared  to  the  experimental  measurements,  the  overall 
pressure  distribution  is  adequately  captured  by  both  fully  turbulent  solutions  over  the  pressure 
side  and  for  the  suction  side  downstream  of  the  region  strongly  influenced  by  laminar-to-turbulent 
transition  at  the  leading  edge  ( x/C  <  0.13).  The  tripless  S-A  results  yield  a  larger  peak  pressure 
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grid 

RANS 

model 

Top  Trip 
(x/C) 

Transition 

(x/C) 

Turbulent 
Separation  (x/C) 

fine  grid 

S-A  tripless 

0.08 

- 

0.88 

DES97 

0,08 

- 

0.75 

DDES 

0.08 

- 

0.86 

coarse  grid  1 

S-A  tripless 

0.08 

- 

6.89 

DES97 

0.08 

- 

0.83 

DDES 

0.08 

- 

0.87 

coarse  grid  2 

DES97 

0.08 

- 

0.82 

DDES 

0.08 

- 

0.87 

Tabic  5:  Comparison  of  S-A  RANS,  DES97  and  DDES  for  a  =  13.3° 


near  the  leading  edge  of  the  suction  side  as  shown  in  figure  17c  as  compared  to  the  fully  turbulent 
runs.  The  larger  peak  suction  is  similar  to  the  behavior  observed  in  the  previous  computations 
which  also  showed  larger  suction  pressures  in  predictions  with  delayed  trailing  edge  separation,  a 
result  attributed  to  changes  in  the  airfoil  circulation  (31),  A  notable  difference  between  the  fully 
turbulent  and  tripless  predictions  is  apparent  for  the  suet  ion -side  pressure  distribution  near  the 
airfoil  leading  edge.  In  the  vicinity  of  x/C  =  0.1,  the  tripless  prediction  results  (cases  with  the  top 
trip  at  x/C  =  0.2  and  x/C  ~  0.13),  show  a  slight  plateau  in  the  pressure  distribution,  a  result  of 
the  laminar  separation  bubble  captured  in  this  calculation,  not  captured  in  the  fully  turbulent  runs. 
Previous  LES  calculations  of  the  flow  were  also  able  to  resolve  the  laminar  separation,  in  addition 
to  details  of  the  transition  process,  though  at  a  substantially  increased  cost  (33). 

The  DES97  and  DDES  pressure  coefficients  are  compared  with  the  S-A  RANS  tripless  results 
in  figure  18.  The  DDES  results  on  different  grids  match  relatively  well  with  the  S-A  RANS  tripless 
results.  While  DES97  results  on  ‘coarse  grid  1  ’  and  ‘coarse  grid  2’  seem  to  be  closer  to  S-A  RANS, 
the  DES97  run  on  the  fine  grid  produces  a  lower  peak  pressure  on  the  suction  side  of  the  airfoil. 

4.2.4  Skin  friction  coefficient 

Figure  1 9  shows  the  skin  friction  coefficient  over  the  suction  side  of  the  airfoil  for  the  RANS  cases. 
The  skin  friction  coefficient  over  the  suction  side  of  the  airfoil  show  that  the  fully  turbulent  cal¬ 
culations  yield  similar  predictions  of  the  skin  friction,  both  S-A  and  SST  in  reasonable  agreement 
with  measurements  for  chordwise  stations  past  x/C  =  0.4.  The  tripless  S-A  results  predict  a  lam¬ 
inar  boundary  layer  prior  to  separation  and  yields  a  substantially  lower  Cf  compared  to  the  fully 
turbulent  runs  upstream  of  x/C  =  0.1.  The  negative  region  in  Cf  near  x/C  =  0.3  identifies  the 
separation  point  and  length  of  the  reverse-flow  region.  For  the  case  in  which  the  top  trip  is  located 
at  x/C  =  0.08,  despite  the  drop  in  Cf,  the  flow  does  not  separate  in  this  region.  Following  reat- 
tachment,  the  skin  friction  in  the  tripless  solution  sharply  increases  to  values  slightly  higher  than  in 
the  fully  turbulent  predictions.  The  trailing  edge  results  show  a  better  agreement  between  the  fully 
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Figure  17:  (a)  Pressure  coefficient  predictions  comparing  tripless  and  fully  turbulent  results;  ( b ) 
Comparison  of  pressure  coefficient  predictions  using  S-A  and  S-A  RC  models;  (c)  Zoomed  view 
of  the  pressure  coefficient  at  the  leading  edge  region;  (rf)  Zoomed  view  of  the  pressure  coefficient 
at  the  trailing  edge  region 
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Figure  18:  Pressure  coefficient:  (a)  Comparison  of  S- A  RANS  and  DES97;  (6)  Comparison  ofS-A 
RANS  and  DDES. 


turbulent  predictions  and  the  measurements  as  compared  to  the  tripless  cases.  Figure  20  shows 
the  comparison  of  DES97  and  DDES  Cf  distributions  with  S-A  RANS  tripless  results.  DES97 
results  show  grid  dependency  with  the  separation  point  moving  further  upstream  in  the  fine  grid 
case  while  DDES  results  are  grid  independent  and  in  good  agreement  with  S-A  RANS. 


Figure  19:  Skin  friction  coefficient  from  the  RANS  cases 


4.2.5  Streamwise  and  wall-normal  velocity  profiles 

The  mean  velocity  components  in  locations  normal  to  the  airfoil  surface  were  measured  in  the  ex¬ 
periments  (30).  Figures  2 1  show  the  comparison  of  the  velocity  profiles  (experiments,  tripless  and 
fully-turbulcnt  RANS)  at  various  stations  on  the  airfoil.  The  predictions  show  the  fully-turbulent  S- 
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Figure  20:  Skin  friction  coefficient:  (a)  Comparison  of  S-A  RANS  and  DES97;  (6)  Comparison  of 
S-A  and  DDES. 


A  and  SST  results  to  be  in  reasonable  agreement  with  the  experimental  measurements  at  x/C  =  0.4 
and  x/C  —  0.5,  with  larger  differences  visible  between  the  experimental  measurements  and  the 
results  closer  to  the  trailing  edge  ( x/C  =  0.775  and  x/C  =  0.87).  The  figures  also  show  that 
the  mean  streamwise  velocity  predicted  in  the  tripless  calculations  lead  that  of  the  fully  turbulent 
solutions. 

Figure  22  show  the  comparison  of  the  streamwise  velocity  profiles  (experiments,  S-A  RANS, 
DES97  and  DDES)  at  various  stations  on  the  airfoil.  The  predictions  show  that  the  DES97  profiles 
near  the  trailing  edge  (at  x/C  =  0.775  and  x/C  =  0.87)  show  discrepancies  when  compared  with 
the  RANS.  Further,  the  DES97  results  on  the  ‘fine  grid’  do  not  match  with  the  results  on  either 
coarse  grid.  In  contrast,  DDES  results  show  no  grid  dependence  and  match  well  with  S-A  RANS. 
Figure  23  show  the  comparison  of  the  streamwise  velocity  profiles  (experiments,  S-A  RANS, 
DES97  and  DDES)  at  various  stations  on  the  airfoil.  These  predictions  indicate  a  similar  trend  with 
the  DES97  profiles  near  the  trailing  edge  (at  x/C  =  0.775  and  x/C  —  0.87)  showing  discrepancies 
when  compared  with  the  RANS,  while  DDES  results  arc  essentially  grid-independent. 

4.2.6  Eddy  viscosity  ratio 

The  eddy  viscosity  ratio  at  two  streamwise  locations  ( x/C  —  0.775  and  x/C  =  0.87)  are  shown 
in  figure  24.  At  x/C  —  0.775,  the  peak  level  of  the  eddy  viscosity  ratio  corresponds  to  about  300 
in  the  RANS  as  well  as  the  DDES  cases.  The  DES97  peak  levels  are  much  lower  with  the  peak 
level  as  low  as  125  in  the  fine  grid  simulation.  This  reduction  in  the  levels  of  the  eddy  viscosity 
ratio  leads  to  MSD.  The  contrast  between  the  DDES  and  DES97  levels  is  even  more  stark  in  the 
eddy  viscosity  profiles  near  the  trailing  edge  of  the  airfoil  at  x/C  -  0.87.  In  the  DDES  cases. 
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Figure  21 :  Streamwise  and  wall-normal  velocity  profiles  at  various  wall-normal  locations  (RANS 
predictions) 


36 


xtt- 0  * 


Ooc» 


S  A  RAM&jVignd) 
DDES  [Pt*  yr«d'| 
DDES  (cont  god  1) 
DDES  [cm™*  gnrf  JJ 


- 


U/U 


0&  09 

u/u 


u/u^  u/u 


Figure  22:  Streamwise  velocity  profiles  at  various  wall-normal  locations  comparing  S-A  RANS, 
DES97  and  DDES 
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Figure  23:  Wall-normal  velocity  profiles  at  various  locations  comparing  S-A  RANS,  DES97  and 
DDES 
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the  boundary  layer  has  been  effectively  shielded  with  RANS  being  maintained  throughout  the 
boundary  layer,  resulting  in  no  reduction  of  eddy  viscosity  and  no  MSD.  This  has  implications  on 
the  separation  point  predicted  by  the  DDES  which  does  not  feel  the  effect  of  modifying  the  grid. 


Figure  24:  Eddy  viscosity  ratio  at  locations  corresponding  to  (a)  x/C  =  0.775;  (5)  x/C  =  0.87. 


4.2.7  Tabulations  of  Lift  and  Drag  coefficients 

The  summary  ofall  the  computations  of  the  airfoil  fora  13.3°  angle  of  attack  are  tabulated  in  tables 
4  and  6.  In  the  experiment,  the  turbulent  separation  location  on  the  suction  side  of  the  airfoil  was 
reported  as  x/C  =  0.83,  and  the  lift  and  drag  coefficients  were  1.56  and  0.0204  respectively.  The 
table  shows  that  the  lift  coefficient  is  higher  in  general  in  the  triplcss  runs  than  in  the  fully  turbulent 
runs.  Among  the  RANS  results  involving  fully  turbulent  runs,  the  S-A  RANS  lift  coefficients 
appear  to  be  closer  to  the  experiments  than  the  SST.  One  of  the  key  points  to  note  is  that  DES97 
produces  a  lift  coefficient  of  1 .40  with  the  fine  grid  and  values  of  1 .52  and  1.51  with  coarse  grids  1 
and  2.  In  contrast  to  the  DES97  results,  the  DDES  lift  coefficients  arc  essentially  grid-independent 
and  also  very  close  to  the  experimental  values. 

4.3  Flow  over  the  Aerospatiale  A-Airfoil  at  a  =  12.3° 

The  RANS  predictions  of  the  two-dimensional  flow  were  computed  using  the  Spalart-Allmaras  (S- 
A)  (with  and  without  rotation  correction)  and  Shear  Stress  Transport  (SST)  turbulence  models.  In 
this  case,  only  the  solutions  of  the  fully-turbulent  flow  were  computed  using  both  RANS  models. 
Even  though  the  pressure  coefficient  in  the  RANS  solutions  are  in  reasonable  agreement  with 
measurements,  the  peak  suction  pressure  is  underpredicted  by  both  models.  The  summary  of  all 
the  computations  of  the  airfoil  for  a  12.3°  angle  of  attack  are  tabulated  in  tables  7  and  the  lift 
and  drag  coefficients  in  table  8.  At  this  angle-of-attack,  all  the  simulations  conducted  were  fully 
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grid 

RANS 

model 

Top  Trip 
location  (i/C) 

cL 

cD 

medium  grid 

S-A  FT 

- 

1.51 

0.032 

S-A  RC  FT 

- 

1.52 

0.032 

SST  FT 

- 

1.43 

0.033 

S-A  tripless 

0.2 

1.62 

0.024 

S-A  RC  tripless 

0.2 

1.62 

0.024 

fine  grid 

S-A  FT 

- 

1.51 

0.033 

S-A  RC  FT 

- 

1.52 

0.032 

SST  FT 

- 

1.45 

0.033 

S-A  tripless 

0.2 

1.62 

0.024 

S-A  RC  tripless 

0.2 

1.62 

0.024 

S-A  tripless 

0.08 

1.60 

0.026 

S-A  tripless 

0.13 

1.62 

0.024 

S-A  tripless 

0.12 

1.62 

0.024 

DES97 

0,08 

1.40 

0.0225 

DDES 

0.08 

1.55 

0.025 

coarse  grid  1 

S-A  tripless 

0.08 

1.60 

0.026 

DES97 

0.08 

1.52 

0.0274 

DDES 

0,08 

1.57 

0.0246 

coarse  grid  2 

DES97 

0.08 

1.51 

0.0226 

DDES 

0.08 

1.57 

0.0248 

Table  6:  Lift  and  Drag  co-efficients  for  the  computations  for  a  =  13.3° 
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turbulenl  RANS.  The  results  obtained  by  the  S-A  RANS,  S-A  RANS  with  rotation  correction  and 
SST  RANS  were  compared  with  the  experimental  results.  The  location  of  turbulent  separation 
predicted  by  the  S-A  RANS  and  S-A  RANS  with  rotation  correction  are  very  close  to  each  other 
(at  x/C  =  0.84  and  xjC  —  0.83  respectively),  while  SST  RANS  predicts  a  turbulent  separation 
point  further  aft,  corresponding  to  x./C  =  0.88.  The  lift  and  drag  coefficients  predicted  are  all  close 
to  each  other  and  compare  reasonably  with  the  experimental  values.  Figure  25  shows  plots  of  the 
pressure  coefficient  distribution  over  the  suction  and  pressure  side  of  the  airfoil.  Figure  26  shows 
the  streamwise  and  wall-normal  velocities  at  points  perpendicular  to  the  wail  at  various  chordwisc 
stations  on  the  airfoil.  The  plots  show  that  all  three  models  predict  the  velocity  profiles  adequately 
at  locations  corresponding  to  x/C  =  0.4  and  x/C  =  0.5,  with  slight  differences  cropping  up 
further  aft,  near  the  location  of  turbulent  separation.  It  appears  that  the  SST  RANS  predicts  the 
wall-normal  velocity  with  greater  accuracy  than  cither  of  the  S-A  RANS  at  x/C  =  0.775  and 
x/C  =  0.87. 


grid 

RANS 

model 

Top  Trip 
location  {x/C) 

Transition 
location  {x/C) 

Turbulenl 
Separation  {x/C) 

Medium  grid 

S-A  FT 

- 

- 

0.84 

S-A  RC  FT 

- 

- 

0.83 

SST  FT 

- 

- 

0.88 

Fine  grid 

S-A  FT 

- 

- 

0.84 

S-A  RC  FT 

“ 

- 

0.83 

SST  FT 

- 

- 

0.88 

Table  7:  Summary  of  the  RANS  computations  for  a  —  12.3° 


grid 

RANS 

model 

Top  Trip 
location  {x/C) 

Cl 

Co 

Medium  Grid 

S-A  FT 

* 

1.42 

0.028 

S-A  RC  FT 

- 

1.44 

0.028 

SST  FT 

- 

1.38 

0.029 

Fine  grid 

S-A  FT 

- 

1.44 

0.029 

S-A  RC  FT 

- 

1.45 

0.028 

SST  FT 

- 

1.39 

0.029 

Table  8:  Lift  and  Drag  co-efficients  for  the  RANS  computations  for  a  =  12.3° 


4.4  Flow  over  the  Aerospatiale  A-Airfoil  at  a  =  7.2° 

The  summary  of  all  the  computations  of  the  airfoil  for  a  7.2°  angle  of  attack  are  tabulated  in  tables 
9  and  the  lift  and  drag  coefficients  in  table  10.  At  thisangle-of-attack,  all  the  simulations  conducted 
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Figure  25:  (n)  Comparison  of  pressure  coefficient  predictions  using  S-A  and  S-A  RC  models; 
(6)  Zoomed  view  of  the  pressure  coefficient  at  the  leading  edge  region;  (<?)  Zoomed  view  of  the 
pressure  coefficient  at  the  trailing  edge  region 
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Figure  26:  Stream  wise  and  wall-normal  velocity  profiles  at  various  wall-normal  locations 
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were  fully  turbulent  RANS.  The  S-A  RANS  and  S-A  RANS  with  rotation  correction  both  predict 
the  same  value  of  the  point  of  turbulent  separation  ( xjC  =  0.97  —  0.98),  while  the  SST  RANS 
fails  to  predict  any  turbulent  separation  on  the  surface  of  the  airfoil.  The  lift  and  drag  coefficients 
predicted  by  all  models  are  close  to  each  other  and  the  experimental  values. 


grid 

RANS 

model 

Top  Trip 
location  (i :jC) 

Transition 
location  (x/C) 

Turbulent 
Separation  (x/C) 

Medium  grid 

S-A  FT 

- 

- 

0.98 

S-A  RC  FT 

- 

- 

0.97 

SST  FT 

- 

- 

- 

Fine  grid 

S-A  FT 

- 

- 

0.97 

S-A  RC  FT 

- 

- 

0.97 

SST  FT 

- 

- 

- 

Table  9:  Summary  of  the  RANS  computations  fora  =  7.2° 


grid 

RANS 

model 

Top  Trip 
location  (j :(C) 

CL 

Co 

Medium  grid 

S-A  FT 

- 

0.98 

0.018 

S-A  RC  FT 

- 

0.98 

0.017 

SST  FT 

- 

0.96 

0.017 

Fine  grid 

S-A  FT 

- 

0.98 

0.017 

S-A  RC  FT 

* 

0.98 

0.017 

SST  FT 

- 

0.96 

0.017 

Table  10:  Lift  and  Drag  co-efficients  for  the  RANS  computations  for  a  =  7.2° 
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Figure  27:  Streamwise  and  wall-normal  velocity  profiles  at  various  wall-normal  locations 
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Figure  28:  (a)  Comparison  of  pressure  coefficient  predictions  using  S-A  and  S-A  RC  models; 
(6)  Zoomed  view  of  the  pressure  coefficient  at  the  leading  edge  region;  (c)  Zoomed  view  of  the 
pressure  coefficient  at  the  trailing  edge  region 
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5  Seeding  boundary  layer  turbulence  in  WMLES 


5.1  Introduction 

Wall-layer  modeling  was  bom  out  of  the  need  to  avoid  the  high  cost  incurred  by  conducting  LES  of 
high  Reynolds  number  turbulent  flows,  wherein,  a  high  computational  cost  results  due  to  the  need 
to  accurately  represent  the  near-wall  eddies.  In  wall-layer  modeling,  the  effects  of  the  eddies  (small 
and  large)  present  in  the  near-wall  region  arc  modeled  in  a  statistical  sense.  This  makes  the  number 
of  grid  points  required  by  the  calculation  of  a  boundary  layer  proportional  to  Re04  (40)  (much 
lower  than  the  grid  requirement  for  a  conventional  LES). 

Wall-layer  models  can  be  generally  classified  into  two  types:  equilibrium  laws  and  zonal  mod¬ 
els  (40).  Equilibrium  laws  assume  the  existence  of  some  generalized  law-of-the-wall.  The  wall 
stress  is  computed  using  this  general  law,  which  is  applied  at  some  distance  from  the  wall.  Even 
though  Equilibrium  taws  have  been  used  with  considerable  success  in  simple  attached  flows,  they 
suffer  from  major  limitations.  For  example,  they  cannot  be  applied  easily  in  complex  geometries 
or  on  fully-unstructured  grids.  A  case  in  point  is  the  calculations  of  the  flow  in  a  rotating  chancel, 
simulated  with  the  use  of  equilibrium  laws  by  Balaras  et  al.  (4).  The  method  adopted  failed  to 
predict  the  quasi-rclaminarization  observed  on  one  side  of  the  channel. 

Zonal  approaches  are  hybrid  RANS/LES  methods  that  solve  the  unsteady  RANS  equations  in 
the  near-wall  region  and  the  filtered  Navier-Stokes  equations  away  from  the  wall.  In  any  zonal 
approach,  there  is  an  explicit  solution  of  a  different  set  of  equations  in  the  inner  layer  and  the 
outer  layer.  For  example,  Two  layer  model  (TLM)  is  a  zonal  approach,  proposed  by  Balaras  and 
Bcnocci  (5)  and  Balaras  et  al.  (4).  The  method  involved  the  solution  of  filtered  Navier-Stokes 
equations  in  the  core  of  the  flow  and  a  simplified  set  of  equations  solved  in  the  wall-layer  (con¬ 
sisting  of  a  fine  grid  embedded  under  a  coarser,  LES  mesh.  The  method  was  subsequently  used  by 
Cabot  (6)  and  Diumo  et  al.  (7)  to  obtain  the  solution  of  the  flow  over  a  backward  facing  step  and 
by  Wang  and  Moin  (8)  in  the  calculations  of  the  trailing  edge  of  the  airfoil.  Conventionally,  DES 
lias  been  used  in  such  a  way  that  the  entire  boundary  layer  is  modeled  by  R  ANS.  Nikitin  et  al.  (72), 
however,  used  DES  as  a  wall-layer  model  in  the  calculations  of  plane  channel  flow.  This  involved 
the  solution  of  the  URANS  equations  in  the  inner  layer  with  the  Reynolds  shear  stress  in  this  region 
being  provided  entirely  by  the  turbulence  model.  In  the  outer  flow,  the  model  was  modified  to  yield 
a  much  lower  eddy  viscosity  as  compared  to  a  pure  RANS  model,  which  allowed  the  formation 
of  turbulent  eddies  capable  of  supporting  most  of  the  Reynolds  shear  stress.  Nikitin  et  al.  (72) 
performed  calculations  with  different  numerical  schemes  and  grids,  over  a  wide  range  of  Reynolds 
numbers  ( Rer  =  180  to  ReT  —  80000).  In  this  study,  it  was  found  that  the  coefficient  of  skin 
friction  was  underpredicted  by  approximately  15%.  The  results  also  featured  the  presence  of  un¬ 
physical  elongated  wall  streaks  in  the  RANS  region  and  at  higher  Reynolds  numbers,  a  logarithmic 
region  with  the  correct  intercept  developed  in  the  quasi-steady  RANS  region  and  an  LES  region 
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characterized  by  an  excessively  high  intercept  of  the  logarithmic  region.  Grid  refinement  only 
moved  the  transition  region  towards  the  wall  without  eliminating  it.  This  discrepancy  between  the 
log-law  intercepts  in  the  inner  and  outer  layer  is  due  to  the  transition  between  the  LES  and  RANS 
regions.  This  region  is  characterized  by  a  lower  eddy  viscosity  while  the  resolved  motions  are  yet 
to  develop.  In  order  to  reach  the  equilibrium  value  of  the  shear  stress,  there  needs  to  be  a  velocity 
gradient  that  balances  the  reduction  in  eddy  viscosity  (40).  Another  method  to  alleviate  this  prob¬ 
lem  involves  the  use  of  a  stochastic  backscatter  model  in  the  inner  layer.  Baggett  (9)  opined  that 
the  stochastic  backscatter  had  the  effect  of  breaking  up  the  elongated  streaks  in  the  RANS  region, 
effectively  de-corrclating  the  velocities  that  enabled  a  rapid  growth  of  resolved  content.  Piomelli  el 
al.  (40)  employed  a  stochastic  backscatter  model  to  improve  the  predictions  of  the  mean  velocity 
profiles  and  skin  friction  coefficient. 

The  use  of  Detached  Eddy  Simulation  (DES)  as  a  wall-layer  model  requires  a  transition  be¬ 
tween  the  RANS  and  LES  regions.  While  this  transition  from  RANS  to  LES  behavior  has  typically 
been  studied  with  evolution  normal  to  the  wall  (4),  (6),  (40).  Of  current  interest  is  the  transition 
from  RANS  to  LES  along  the  main  flow  direction,  which  can  result  in  applications  because  of 
streamwise  grid  refinement.  In  either  scenario,  a  DES  solution  is  characterized  by  a  reduction 
in  the  modeled  length  scale  which  draws  down  the  eddy  viscosity,  thereby  lowering  the  modeled 
stress.  This  reduction  in  modeled  stress  requires  a  corresponding  increase  in  the  resolved  stress 
near  the  transition  layer  in  order  to  avoid  degrading  predictions  of  the  skin  friction  and  the  mean 
velocity  profile. 

The  early  calculations  of  temporally  developing  flows  (homogeneous  isotropic  decay,  shear 
layers,  mixing  layers,  and  plane  channels)  involved  the  use  of  periodic  boundary  conditions  in 
the  flow  direction.  These  conditions  supply  a  realistic  turbulent  field  at  the  inflow,  where  physi¬ 
cally  realistic  turbulence  can  be  introduced  by  recycling  the  outflow  plane.  In  spatially  developing 
flows,  on  the  other  hand,  inflow  conditions  must  be  assigned  in  a  manner  consistent  with  the  basic 
requirements  and  properties  of  LES.  The  early  efforts  in  specifying  inflow  boundary  conditions 
involved  the  modification  of  periodic  conditions  to  supply  the  inflow  velocity.  Spalart  (45)  per¬ 
formed  DNS  of  sink-flow  boundary  layers  and  flat-plate  boundary  layers  using  periodic  boundary 
conditions.  To  take  into  account  the  variable  boundary  layer  thickness,  source  terms  were  used  in 
the  equations  to  transform  them  into  a  selfsimilar  coordinate  flow  in  which  the  flow  was  periodic. 
Spalart  and  Watmuff  (10)  modified  this  method  to  include  a  "fringe”  region  appended  to  the  end  of 
the  domain,  in  which  forcing  terms  were  added  to  the  momentum  equations  to  decrease  the  bound¬ 
ary  layer  thickness  and  re-establish  an  equilibrium  boundary  layer.  Periodic  boundary  conditions 
were  then  used  to  reintroduce  the  outlet  velocity  field  at  the  inlet.  Lund  et  ai  citclund  proposed 
a  technique  developed  for  flat  plate  boundary  layers,  in  which  a  plane  of  data  from  a  location  sev¬ 
eral  boundary-layer  thicknesses  d  downstream  of  the  inflow  plane  was  manipulated  by  rescaling 
the  inner  and  outer  layers  of  velocity  profiles  separately,  to  account  for  the  different  scaling  laws 
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observed  in  these  regions.  The  rescaled  profiles  were  then  re-introduced  at  the  inlet. 

Inflow  conditions  can  also  be  generated  by  running  a  separate  precursor  calculation  of  an  equi¬ 
librium  flow  in  which  periodic  boundary  conditions  can  be  used.  The  velocity  field  in  one  plane 
normal  to  the  streamwise  direction  can  be  stored  at  each  time-step.  The  sequence  of  planes  is  then 
read  in  as  inflow  data  for  a  separate  calculation.  This  method  has  been  used  with  considerable 
success  in  the  studies  of  Li  et  ai  (11)  and  Piomelli  et  al.  ( 1 2).  Another  class  of  methods  that  do 
not  involve  recycling  are  generally  based  on  the  generation  of  synthetic  turbulence  with  assigned 
moments  and  spectra,  generated  using  the  manipulation  of  random  sequences.  Le  et  ai  ( 1 3)  per¬ 
formed  calculations  of  a  backward  facing  step  in  which  a  mean  velocity  profile  was  assigned  at  the 
inflow  plane  to  which  random  fluctuations  were  superposed.  To  determine  the  random  fluctuations, 
they  used  the  method  proposed  by  Lee  et  ai  (14).  In  this  study,  it  was  found  that  the  turbulence 
levels  decayed  rapidly  (13)  because  of  the  lack  of  phase  information  of  the  real  turbulent  eddies. 

Batten  et  ai  (69)  proposed  a  new  method  to  generate  synthetic  turbulence,  that  takes  into  ac¬ 
count  the  anisotropy  of  the  flow.  The  method  is  based  on  the  superposition  of  sinusoidal  modes 
with  random  frequencies  and  wavenumbers,  with  given  moments  and  spectra.  The  approach  in¬ 
volved  the  modification  of  wavenumbers  to  yield  eddies  that  are  more  elongated  in  the  direction  of 
larger  Reynolds  stresses,  thereby  introducing  more  realistic  anisotropic  eddies  into  the  flowfield. 
A  summary  of  the  synthetic  turbulence  method  is  provided  below:  These  investigators  applied 
their  method  to  a  hybrid  RANS/LES  simulation  of  turbulent  channel  flow  and  found  that  resolved 
Reynolds  shear  stress  decreased  downstream  of  the  inflow.  While  the  shear  stress  levels  sub¬ 
sequently  increased  with  continued  downstream  evolution,  the  expected  levels  were  lower  than 
anticipated  even  ten  channel  half-heights  downstream  of  the  inlet  boundary.  This  discrepancy  was 
attributed  to  an  incorrect  behavior  of  the  pressure-strain  term  in  the  Reynolds  stress  budget  (69), 
responsible  for  rapid  decay  of  the  wall-normal  ((uV))  velocity  variance,  resulting  in  lower  values 
of  resolved  Reynolds  shear  stresses.  To  more  quickly  re-establish  the  correct  Reynolds  shear  stress 
levels,  Spille-KohofT  and  Kaltcnbach  (50)  proposed  a  method  in  which  a  synthetic  turbulent  field 
is  prescribed  at  the  inflow  plane  and  a  number  of  control  planes  are  employed  downstream  of  the 
inlet.  These  control  planes  amplify  the  wall-normal  velocity  fluctuations  in  order  to  drive  the  re¬ 
solved  shear  stress  towards  a  target  profile.  Piomelli  et  al.  (47)  applied  this  scheme  to  simulations 
of  a  spatially-developing  channel  flow. 

In  this  study  (47),  Proportional-Integral  (PI)  controllers  were  used  to  increase  the  resolved 
stress  levels  towards  the  target  profile.  A  PI  controller  is  a  generic  control  loop  feedback  mecha¬ 
nism  widely  used  in  industrial  control  systems.  It  attempts  to  correct  the  error  between  a  measured 
process  variable  and  a  desired  setpoint  (desired  output),  by  calculating  and  then  initiating  a  cor¬ 
rective  action  that  can  adjust  the  process  accordingly.  The  PI  controller  algorithm  involves  the 
interplay  between  two  values:  the  proportional  values  and  the  integral  values.  The  proportional 
value  determines  the  reaction  to  the  current  error  while  the  integral  value  determines  the  rcac- 
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lion  based  on  the  sum  of  recent  errors.  The  weighted  sum  of  these  two  actions  is  ucd  to  adjust 
the  controlled  variable.  PI  controllers  with  proportional  and  integral  gains  corresponding  to  1 
and  50,  respectively,  were  used  to  amplify  the  wall-normal  fluctuations  downstream  of  the  inflow 
plane  (47).  The  study  (47)  found  that  the  error  in  the  Reynolds  shear  stress  was  reduced  to  ac¬ 
ceptable  levels  within  five  channel  half-heights  while  the  skin  friction  coefficient  required  longer 
recovery  lengths  (around  15  channel  half-heights  downstream  of  the  inflow  plane)  (47).  Keating  et 
al.  (61)  incorporated  a  similar  synthetic  turbulence  method  to  supply  inflow  boundary  data,  but 
with  reduced  information  in  the  assignment  of  moments,  and  assessed  the  effect  on  the  flow  field 
predictions.  These  investigators  modified  the  forcing  method,  by  gradually  increasing  the  force 
over  a  number  of  planes  upstream  of  each  control  plane,  ensuring  that  the  resolved  Reynolds  shear 
stress  reached  target  levels  using  two  control  planes,  spaced  !  .3  boundary  layer  thicknesses  apart. 
The  transient  period  for  the  controllers  was  significantly  reduced  (without  adverse  effects  on  other 
parameters)  by  using  values  of  proportional  and  integral  gains  corresponding  to  1  and  30  respec¬ 
tively  (61).  A  description  of  the  synthetic  turbulence  method  is  provided  in  the  following  section. 


5.2  Synthetic  turbulence 


The  use  of  synthetic  turbulence  requires  the  specification  of  complete  statistical  data  (full  Rcynolds- 
slress  tensor  and  dissipation  rate).  Typically  when  the  S-A  model  is  used  as  the  RANS  model,  these 
details  arc  not  available  and  this  section  describes  the  use  of  the  synthetic  turbulence  method  to 
generate  the  inflow  condition  in  the  presence  of  reduced  information  (61).  The  method  consists  of 
the  manipulation  of  sines  and  cosines  with  random  phases  and  amplitudes  to  obtain  an  intermediate 
velocity  field: 


p*cos(djXj  +  unt)  +  q"sin(djX}  4-  ujnt) 


(40) 


where, 

Xj  =  27TXj/Lb  (41) 

t  =  27t  t/rb  (42) 

are  spatial  coordinates  normalized  by  the  length-  and  timescale  of  turbulence;  rj  —  k.j<  and  Lb  — 
TbVb  are  the  turbulence  timescale  and  turbulence  lengthscale,  and  Vj,  is  given  by  k where  k  is 
the  turbulent  kinetic  energy.  The  random  frequencies  wn  =  N( 1. 1)  are  taken  from  a  normal 
distribution  N(fi,  rr2)  with  a  mean  and  variance  of  1.  The  amplitudes  are  given  by. 


p”  —  CijfcCj  dk  (43) 

9,"  =  iijkZjndkn  (44) 
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where  Cj"  and  are  given  by  N(0. 1),  and. 


(45) 


are  modified  wavenumbers  obtained  by  multiplying  the  wavenumbers  d,n  by  Vb/cn,  where  cn  is 
given  by, 


cn  = 


_  /u  <u  ') 

2  '  m  '  d) Jt/J? 


(46) 


The  isotropic  distribution  of  wavenumbers  d3  is  thus  transformed  into  an  anisotropic  one,  result¬ 
ing  in  eddies  stretched  in  the  directions  of  greatest  Reynolds  stress.  The  wave-numbers  d,n  = 
N( 0, 1/2)  are  chosen  from  a  normal  distribution  with  variance  1/2.  A  tensor  scaling  is  used  to 
reconstruct  the  synthetic  turbulent  fluctuation  field, 


Ui  — 


(47) 


where  a**  is  the  Chotesky  decomposition  of  the  Reynolds  stress  tensor.  When  only  reduced  in¬ 
formation  is  available,  the  method  can  be  implemented  as  described  below.  When  the  S-A  RANS 
turbulence  model  is  used,  the  eddy  viscosity  is  related  to  the  turbulent  kinetic  energy  using  the 
experimental  result  of  Bradshaw  et  al,  (1 ). 


|-  (uV)|  =  t/, 


du 

Oy 


=  aik 


(48) 


where  a(  —  y/c^  and  cM  —  0.09.  To  obtain  the  individual  stresses  ( u'u '),  {v'v')  and  {w'w'}, 
the  ratio  of  the  above  three  stresses  for  a  nominal  boundary  layer  is  obtained  (from  experimental 
observations  or  Direct  Numerical  Simulations  (DNS)).  In  the  current  study,  the  ratio  of  the  three 
stresses  were  obtained  from  DNS  of  a  channel  flow  al  ReT  =  180  (3).  Once  these  ratios  arc 
obtained,  the  levels  of  the  individual  stresses  that  add  up  to  form  k,  calculated  using  48  can  be 
obtained.  To  obtain  the  timescale,  the  definition  of  eddy  viscosity  from  the  k  —  (  two-equation 
turbulence  model  can  be  used, 

fc2 

t  =  c»~  (49) 

k  1 
n  =  -  =  — 


HU 

dy 


5.3  Boundary  layer  stirring  and  control  methodology 

As  discussed  earlier,  the  use  of  DES  as  a  wall-layer  model  results  in  a  transition  zone  between  the 
RANS  and  the  LES  regions,  characterized  by  a  lower  modeled  stress  and  requiring  an  increase  in 
the  resolved  stress  near  the  transition  layer.  The  resolved  Reynolds  stress  could  be  increased  by 
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providing  realistic  fluctuations  in  the  transition  layer,  by  seeding  the  boundary  layer  with  three- 
dimensional  content,  to  achieve  a  result  similar  to  that  accomplished  by  the  use  of  inflow  boundary 
conditions.  In  the  current  effort,  resolved  stresses  are  generated  by  seeding  the  near- wall  flow  with 
fluctuations  in  order  to  generate  three-dimensional  structure  in  the  boundary  layer.  Downstream 
of  the  streamwise  location  at  which  the  fluctuations  are  seeded,  control  planes  are  employed  to 
amplify  the  wall-normal  velocity  fluctuations  in  order  to  tune  the  resolved  shear  stress  levels  based 
on  a  target  profile.  In  this  part  of  the  study,  a  stochastic  force  added  to  the  momentum  equations 
is  used  to  generate  velocity  fluctuations.  These  fluctuations  are  subsequently  amplified  based  on 
a  prescribed  shear  stress  profile  at  a  number  of  control  planes  downstream  of  the  stochastically 
stirred  region.  Velocity  fluctuations  were  seeded  into  the  near-wall  flow  using  a  stochastic  force 
formed  from  Gaussian  random  numbers  scaled  such  that  the  magnitude  of  the  force  is  maximum 
near  the  location  of  the  RANS-LES  interface.  Equation  (5 1 )  shows  one  of  the  forms  of  the  stirring 
force  used  to  generate  fluctuations: 


fs,i  — 


l  +  (Ay)27s'’‘ 


(51) 


In  (51),  A  is  a  constant  that  ensures  that  the  amplitude  of  the  stirring  force  is  maximum  near  the 
RANS-LES  interface;  As  represents  the  amplitude  of  the  stirring  force.  Preliminary  investigations 
were  also  conducted  using  different  forms  for  the  stirring  force  based  on  the  values  of  Aa  (in  each 
of  the  three  directions)  and  A  in  order  to  assess  the  influence  on  the  resulting  velocity  fluctuations. 


Inlet  stirring  control  Developed  turbulent 

section  region  planes  boundary  layer 


Figure  29:  Schematic  depicting  the  inlet,  stirring  region,  and  controllers.  The  stirring  region  is 
applied  over  a  streamwise  distance  of  0.55  and  four  control  planes  are  spaced  5  apart.  The  distance 
between  the  end  of  the  stirring  region  and  the  first  control  plane  is  0.25. 

The  controller  uses  the  error  in  the  Reynolds  shear  stress  as  its  input  parameter.  The  error  is 
obtained  by  a  comparison  of  the  Reynolds  shear  stress  to  a  prescribed  set  point  (47): 

e{y,  t )  =  (uV)  (x0,  y )  -  (uV)*,(  (x0,  V •  t)  (52) 

where  (uV)  (xq.  y)  is  the  target  Reynolds  stress  at  streamwise  reference  location  x  =  x0,  and 
(u'v'y  1  (xq.  y.  t)  is  the  current  value,  averaged  over  the  spanwise  direction  and  over  time  using  a 
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short  backward  time  average, 


{u'v'Y'1  [t  +  A t)  =  ^~  {u'v'Y  4-  ( 1  -  (uV)*'1  ( t ) .  (53) 

In  (53),  Tavg  ~  25 /Uc,  the  superscripts  2  and  f  indicate  averaging  over  the  spanwise  direction  and 
time,  and  Uc  refers  to  the  velocity  at  the  center  of  the  channel.  As  described  in  Piomelli  et  al.  (47) 
and  Kaltcnbach  et  at.  (50),  the  forcing  aims  at  the  enhancement  or  damping  of  local  events  that 
contribute  to  the  Reynolds  shear  stress.  This  is  achieved  by  the  formula, 

Ic(x 0,  y,  2,  t)  =  r(y,  t)  (u(x0,  V,  M)  -  {uY‘  (x0 ,y))  ,  (54) 

where  the  magnitude  is  controlled  by  the  error  according  to, 

r(y,  t)  ~  Kpe{y ,  t)  +  K,  f  e(y,  a)da ,  (55) 

JO 

where  Kp  and  Kj  are  the  proportional  and  integral  gains  of  the  controllers,  respectively.  Pi¬ 
omelli  et  al.  (47)  reported  that  unrealistically  high  shear  stress  events  were  prevented  by  applying 
the  control  force  only  when  |u'|  <  7 U14  and  |u'|  <  7 „(/*,  where  7„  =  0.6,  7„  =  0.4,  and  Ub 
is  the  bulk  velocity.  Further,  forcing  employed  (47)  was  concentrated  on  more  energetic  events 
by  applying  the  force  only  when  the  additional  criterion  |uV|  >  0.00016(7^  was  satisfied  (47). 
These  additional  conditions  were  not  used  in  the  control  schemes  implemented  in  the  current  work 
as  discussed  below.  The  most  significant  issue  in  the  application  of  the  method  is  that  the  use 
of  controllers  (i.c.,  body  forces)  requires  the  specification  of  gain  values  that  dictate  the  rate  of 
convergence  of  the  resolved  shear  stress  to  the  target  values  and  also  ensure  that  the  forcing  ap¬ 
plied  is  not  too  large  or  unphysical.  If  chosen  incorrectly,  the  gain  values  induce  instabilities  in  the 
control  process,  resulting  in  an  output  that  diverges.  Even  though  this  can  be  rescinded  by  the  use 
of  saturations  which  deactivate  the  control  commands  and  ensure  that  the  forcing  is  not  applied 
whenever  a  predefined  criterion  is  not  satisfied,  this  is  not  desirable  as  it  has  implications  on  the 
control  process.  The  primary  focus  of  the  current  effort  is  to  investigate  the  control  process  and  to 
analyse  the  effect  of  different  control  methodologies  (use  of  Proportional- Integral  controllers,  use 
of  integral  controllers,  use  of  non-linear  Proportional-Integral  controllers  based  on  non-linearities 
in  the  system  and  adaptive  controllers)  in  achieving  the  desired  levels  of  resolved  shear  stress  at 
the  specific  control  planes. 

As  described  below,  early  stages  of  the  investigation  involved  the  system  identification  process, 
with  the  objective  of  identifying  a  linear  model  formed  from  the  open  loop  response  of  the  resolved 
shear  stress  using  random  inputs  at  the  control  planes.  Optimal  gains  for  the  controller  are  spec¬ 
ified  based  on  the  location  of  the  poles  of  the  linear  model  in  the  closed-loop  configuration.  The 
shortcomings  of  this  method  are  explained  in  detail  later  in  this  chapter. 

During  the  next  stage  of  the  investigation,  the  response  of  the  non-linear  system  to  preset  gain 
values  was  measured  during  a  precursor  simulation.  These  gain  values  were  then  incremented 
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unlit  the  onset  of  numerical  instabilities.  The  value  of  gain  obtained  prior  to  the  onset  of  numerical 
instabilities  can  then  be  used  in  the  final  simulations.  It  is  to  be  noted  that  such  a  method  is 
used  to  arrive  at  the  gain  for  an  integral  controller  and  has  been  described  in  greater  detail  in 
the  later  sections.  This  method  was  not  found  to  work  accurately  despite  the  use  of  very  small 
increments  of  the  gain  values  {of  the  order  of  10~6)  in  the  precursor  simulation.  It  was  found  that 
the  during  the  initial  stages  of  the  simulation,  due  to  the  absence  of  any  fluctuations  the  gain  values 
would  increase  and  reach  higher  values  comprimising  the  stability  of  the  precursor  simulation. 
These  values  of  gains  were  not  small  enough  to  ensure  the  stability  of  the  simulation  of  interest. 
Subsequently,  a  deep  investigation  of  the  control  process  and  the  effect  of  the  use  of  different 
controllers  on  the  system  was  conducted.  In  this  part  of  the  work,  the  behavior  of  the  error  in  the 
system  for  various  controllers  was  thoroughly  investigated.  It  was  found  (hat  the  wall  shear  stress 
and  resolved  Reynolds  stress  levels  are  independent  of  the  type  of  controller  used  in  the  method, 
implying  that  there  is  little  to  be  gained  in  the  use  of  a  PI  controller  over  an  integral  controller. 
However,  the  specification  of  optimal  gains  that  keep  the  system  stable  is  a  key  issue  and  to  make 
the  process  less  ad-hoc,  an  adaptive  integral  controller  was  developed.  This  adaptive  controller 
adjusts  the  integral  gain  value  during  the  simulation  ensuring  that  the  stability  of  the  system  is  not 
compromised. 

The  use  of  integrator-based  controllers  to  achieve  target  resolved  stress  levels  has  also  been 
tested  in  channel  flow  with  the  use  of  the  ‘synthetic  turbulence’  inflow  condition  of  Batten  el  al 
(69).  In  this  case,  the  three-dimensional  content  is  supplied  by  the  inflow  condition,  and  the  use  of 
the  stochastic  force  in  the  momentum  equations  is  obviated.  Some  of  these  predictions  have  also 
been  included  in  this  report.  Other  issues  involved  in  the  method,  such  as  the  effect  of  the  form 
and  magnitude  of  the  stirring  force  on  the  mis  values  of  the  velocity  fluctuations,  the  efTect  of  grid 
refinement,  and  number  of  controllers  used  on  the  wall  shear  stress  are  also  discussed. 

5.4  Simulation  overview 

5.4.1  Numerical  method 

Spatial  derivatives  in  the  governing  equations  were  discretized  using  second-order  central  differ¬ 
ences  on  a  staggered  grid.  In  the  momentum  equations,  the  diffusion  terms  are  advanced  using 
Crank-Nicolson,  while  the  remaining  terms  are  advanced  explicitly  using  a  second-order  Adams- 
Bash forth  method.  For  the  transport  equation  for  the  variable  v  in  the  S-A  mode!  equation,  the 
diffusion  term  is  treated  implicitly  while  all  other  terms  are  treated  explicitly.  The  equations  are 
advanced  in  time  using  a  semi-implicit  fractional  step  method.  The  Poisson  equation  for  pressure 
is  solved  using  a  direct  approach  based  on  Fourier  transforms.  No-slip  boundary  conditions  are 
applied  at  the  lower  and  upper  walls  of  the  channel.  Periodic  boundary  conditions  were  used  in 
the  spanwise  direction  and  the  inlet  boundary  condition  was  fixed  with  a  RANS  profile  making-up 
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the  specified  variables.  A  convective  boundary  condition  was  employed  at  the  outlet  of  the  domain 
with  the  outflow  velocity  corrected  at  every  time  step  to  ensure  mass  conservation. 

5.4.2  Computational  grid 

Computations  were  performed  at  Reynolds  numbers  of  400  and  5000,  based  on  friction  velocity 
and  channel  half-width  S.  The  grid  sizes  used  for  all  the  simulations  are  summarized  in  Table  1 1 . 
The  spanwise  extent  of  the  channel  was  fixed  at  n  for  all  the  cases.  For  the  coarse  grid  with 
=  0,1,  the  RANS-LES  interface  was  located  at  a  wall-normal  distance  y/S  =  0.065  from  the 
channel  walls  (corresponding  to  y+  =  26  when  Rer  =  400  and  y*  =  325  when  ReT  —  5000). 
For  the  fine-grid  runs  with  the  grid  size  129  x  75  x  65,  the  number  of  controllers  and  the  spacing 
between  the  controllers  was  also  varied  to  study  their  effect  on  the  wall  shear  stress.  The  runs 
conducted  have  been  summarized  in  Table  12. 


Rer 

grid 

Method 

A, 

A, 

Streamwise  extent 

400 

129  x  75  x  33 

Boundary  layer  stirring 
with  controlled  forcing 

0.1 

0.1 

4tt 

129  x  75  x  33 

Synthetic  turbulence  with 
controlled  forcing 

0.1 

0.1 

4tt 

259  x  75  x  33 

Synthetic  turbulence  with 
controlled  forcing 

0.1 

0.1 

8?r 

5000 

129  x  75  x  33 

Boundary  layer  stirring 
with  controlled  forcing 

0.1 

0.1 

4tt 

129  x  75  x  65 

Boundary  layer  stirring 
with  controlled  forcing 

0.1 

0.05 

Ait 

195  x  75  x  65 

Boundary  layer  stirring 
with  controlled  forcing 

0.1 

0.05 

Cm 

259  x  75  x  65 

Boundary  layer  stirring 
with  controlled  forcing 

0.1 

0,05 

8tt 

Table  1 1 :  Grid  sizes  used  for  the  channel  flow  investigations. 


Number  of  controllers 

Spacing  between  controllers 

5 

x/S  -  1.0 

5 

x/6  =  0.5 

10 

t/6  =  0.5 

Table  12:  Runs  conducted  with  the  fine  grid  with  the  grid  size  129  x  75x33  cells  in  the  streamwise, 
wall-normal,  and  spanwise  directions,  respectively. 
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5.5  Design  of  the  controllers 

PI  controllers  have  widespread  applications  in  process  control,  particularly  when  the  mathematical 
model  of  the  system  is  unknown  and  analytical  methods  for  system  design  are  not  applicable.  The 
PI  controller  shown  in  (55),  applied  by  Piomelli  et  ai  (47)  has  the  standard  form,  consisting  of 
a  proportional  term  and  an  integral  term.  The  proportional  term  ensures  that  the  controller  output 
(r(y,  (})  is  proportional  to  the  error  ( e(y ,  f))  of  the  system  while  the  integral  term  ensures  that  the 
error  accumulated  over  a  period  of  time  is  taken  into  account  by  multiplying  this  quantity  to  a 
constant  and  then  added  to  the  controlled  quantity  (r{y,  i)). 

The  open  loop  transfer  function  (ratio  of  the  2-transform  of  the  output  to  the  2-transform  of 
the  input)  of  the  PI  controller  in  the  discretized  form  is  given  by  (56),  where  Kp  represents  the 
proportional  gain  and  Kt  represents  the  integral  gain. 

C(j)  =  Kp  I  1  +  T P— ^ - -  1  (56) 

The  PI  controller  described  above  forms  a  part  of  a  non-linear  system,  the  solution  of  which  can 
be  considered  to  be  the  whole-integration  scheme  of  the  Navier-Stokes  equations  (inclusive  of  the 
body  forces). 

A  detailed  analysis  of  the  control  process  has  been  carried  out  with  the  goal  of  maximizing  the 
rate  of  convergence  of  the  resolved  stress  to  the  desired  values  and  also  ensure  the  stability  of  the 
closed  loop  system  without  the  use  of  any  ad-hoc  adjustments  in  order  to  prevent  the  numerical 
method  from  diverging.  'Stability’  in  this  case  implies  that  the  output  shear  stress  remains  bounded 
over  any  amount  of  time  for  a  bounded  input  control  force. 

Prior  to  the  discussion  on  system  identification  of  a  linear  model,  results  showing  that  the 
a  purely  integral  controller  is  sufficient  in  achieving  the  desired  results  arc  shown.  Figure  30 
shows  the  responses  of  the  system,  i.e.,  time  development  of  (u'v')2,1,  at  the  third  control  plane. 
Shown  arc  time  histories  at  two  wall-normal  locations  between  the  wall  and  RANS-LES  interface, 
y/S  =  0.0321  ( y+  —  161)  and  y/5  =  0.0168  ( y+  —  84)  along  with  the  set  point  for  a  case  with 
Kp  —  0  and  A"/  =  0.5.  In  the  figures,  filtered  response  refers  to  a  low  pass  filtered  form  of 
the  output  that  confirms  that  the  mean  of  the  output  is  very  close  to  the  set  point.  The  mean  of 
the  response  matches  the  set  point  very  closely,  showing  that  the  integral  controller  is  effective  in 
eliminating  the  steady  slate  error. 

5.5.1  Identification  of  a  linear  model 

Initial  efforts  were  focused  on  the  identification  of  a  linear  model,  and  arriving  at  integrator  gain 
values  based  on  the  model  poles.  During  the  identification  step,  controller  parameters  arc  chosen 
so  as  to  produce  outputs  (i.e.,  resolved  shear  stress)  that  are  of  the  same  order  of  magnitude  as 
the  desired  outputs;  this  ensures  that  the  system  operates  in  the  vicinity  of  the  set  point  (the  target 
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Figure  30:  System  response  at  wall-normal  locations  yft 5  —  0.0321  (y+  =  161)  and  y/ri  =  0.0168 
(y+  =  84)  at  the  third  control  plane  x/6  —  2.15. 

Reynolds  shear  stress).  In  this  first  step,  Gaussian-distributed  noise  (x)  is  added  by  each  controller 
instead  of  the  control  force,  and  the  inputs  and  corresponding  outputs  (resolved  shear  stress)  are 
recorded.  A  linear  model  is  fit  to  these  data,  and  since  the  data  have  been  collected  in  the  vicinity 
of  the  desired  set  point  for  the  system,  the  resulting  model  can  be  considered  as  an  approximation 
of  the  linearization  of  the  numerical  scheme  around  the  set  point.  With  such  a  setup,  the  system  is 
said  to  be  in  open-loop  mode  as  there  is  no  control  action  that  determines  the  inputs  that  are  added. 
The  state  space  representation  of  the  linear  model  is  given  below: 

xe+i  =  Axt  4-  B\  (57) 

where, 

•  x(  represents  the  difference  between  the  output  of  the  system  -  (uV)1,  (f )  and  the  setpoint; 

•  x  represents  the  input  to  the  system  -  the  noise  added  in  the  open  loop  system; 

•  A  and  B  are  the  state  matrix  and  input  matrix  respectively,  with  A  representing  the  slate  of 
the  system  (matrix  that  represents  the  dynamics  of  the  system,  that  summarizes  the  effect  of 
past  inputs  on  future  outputs  of  the  system). 

Studying  the  dynamics  of  the  system  in  the  open-loop  configuration  give  insights  on  the  underlying 
structure  of  the  system,  accomplished  by  examining  the  eigenvalues  of  the  dynamic  matrix  A. 
These  eigenvalues  represent  the  poles  of  the  system  under  consideration.  Equation  57  can  be 
manipulated  to  obtain  the  matrices  ,4  and  B, 

X  =  \A  B\\X,  x\T 
X[X,  x ]  =  [/l  B\\X,  x)T\X.  x\ 

\A  B)  =  X[X,  *]([*,  x\T\X,  ,|) 


where, 
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(58) 

(59) 

(60) 


Figure  3 1 :  eigenvalues  of  the  open-loop  configuration  (a)  with  rj  approximately  700%  of  the  output 
of  the  controller;  (b)  with  77  approximately  10%  of  the  output  of  the  controller. 

*  X  represents  a  vector  that  contains  the  difference  between  the  output  of  the  system  -  (u'v') 
and  the  set  point  at  time  level  n  +  1  (from  data  collected  during  the  identification  step); 

•  Xs  represents  a  vector  that  contains  the  difference  between  the  output  of  the  system  -  {u'v') 
and  the  set  point  at  time  level  n  (from  data  collected  during  the  identification  step). 

Figures  31a,b  show  the  eigenvalues  of  the  dynamic  matrix  (A)  for  two  cases  in  the  open-loop 
configuration.  Figure  31a  shows  the  eigenvalues  of  the  dynamic  matrix  for  a  case  in  which  the 
standard  deviation  of  the  noise  r;  was  approximately  700%  of  the  output  resolved  shear  stress.  Such 
a  large  magnitude  of  the  standard  deviation  of  the  noise  rj  enables  a  more  robust  characterization. 
Figure  3 1  b  shows  the  eigenvalues  of  the  dynamic  matrix  for  a  case  in  which  the  standard  deviation 
of  the  noise  \  was  approximately  10%  of  the  output  of  the  controller.  The  figures  show  that  the 
eigenvalues  of  the  dynamic  matrix  A  are  invariably  close  to  the  edges  of  the  unit  circle.  Thus, 
system  identification  can  be  accomplished  using  x  that  has  a  much  higher  standard  deviation  than 
the  term  r(y,  f)  in  equation  55,  increasing  the  robustness  of  the  characterization. 

5.5.2  Closing  the  loop  and  obtaining  the  gain  values 

A  system  is  said  to  be  in  a  closed-loop  configuration  when  there  is  an  active  feedback  control  loop 
that  acts  to  control  the  output  of  the  system.  The  closed  loop  response  of  the  system  depends  upon 
the  location  of  the  poles  of  the  linear  model  in  a  closed  loop  configuration.  Poles  and  zeros  of 
a  transfer  function  are  the  frequencies  for  which  the  transfer  function  becomes  infinity  and  zero 
respectively.  The  values  of  the  poles  and  zeros  of  a  system  determine  if  the  system  is  stable  and 
the  response  of  the  system.  Control  systems  can  be  designed  by  assigning  specific  values  to  the 
poles  and  the  zeros  of  the  system.  The  loop  can  be  closed  by  the  placement  of  a  pole  at  the  edge 
of  the  unit  circle  corresponding  to  z  —  1  and  arriving  at  an  appropriate  gain  value  for  the  integral 
controller.  Since  it  is  desired  to  have  the  system  output  match  the  set  point  exactly,  and  since 
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Figure  32:  Root  Locus  Plot  for  a  particular  combination  of  output  and  input.  In  this  diagram,  o 
represents  the  zeros  of  the  system  and  x  represents  the  poles  of  the  system. 

the  set  point  is  constant  with  time,  then  the  integral  controller  itself  is  sufficient  for  this  purpose. 
The  proportional  term  is  therefore  not  useful,  thus  tl  can  be  set  to  zero.  The  results  shown  below 
confirmed  this  conclusion,  by  comparing  the  effect  of  a  controller  using  the  integral  term  only  with 
other  controllers  using  different  values  of  Kp.  An  attempt  was  made  to  find  the  optimal  value  of  the 
integrator  gain  by  studying  the  locus  of  the  poles  of  the  transfer  function  of  the  system  in  closed* 
loop  as  the  integrator  gain  is  varied.  To  arrive  at  the  optimal  gain,  the  MIMO  (Multiple-Input 
Multiple-Output)  model  was  split  into  a  number  of  S1SO  (Single-Input  Single-Output)  systems 
and  the  smallest  value  of  the  gain  for  which  the  model  remains  stable  was  selected.  The  optimal 
gain  is  the  minimum  value  among  the  gains  arrived  at  for  all  the  sets  of  individual  inputs  and 
outputs  that  ensures  that  the  system  remains  stable. 

The  threshold  of  stability  for  the  system  is  defined  on  the  basis  of  the  location  of  the  poles  of 
the  model.  The  definition  of  the  threshold  is  justified  because  the  corresponding  system  is  stable 
in  the  open-loop  configuration  despite  the  model  having  poles  outside  the  unit  circle.  This  is  not 
a  contradiction,  because  the  linear  model  approximates  the  behavior  of  an  underlying  nonlinear 
system,  which  is  the  numerical  integration  scheme.  Therefore,  we  can  consider  the  circular  bound¬ 
ary  marked  by  the  outermost  linear  model  poles  as  a  “safety  region",  in  the  sense  that  we  know 
that  the  corresponding  nonlinear  system  is  stable  for  that  given  location  of  the  outermost  linear 
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Figure  33:  Variation  in  the  gain  values  as  the  threshold  of  stability  is  varied. 

model  poles.  We  can  also  assume  that  if  the  poles  of  the  controlled  linear  model  remain  within 
that  region,  then  the  corresponding  controlled  nonlinear  system  will  also  be  stable.  On  the  other 
hand,  since  the  poles  of  the  controlled  model  are  farther  from  the  origin  and  closer  to  the  edges  of 
the  unit  circle,  the  controlled  system  will  have  a  more  rapid  response,  i.e.,  faster  convergence  of 
the  resolved  shear  stress  towards  the  desired  values.  Figure  32  shows  a  representative  plot  of  the 
root  locus  for  a  particular  combination  of  output  and  input.  The  plot  also  shows  the  threshold  of 
stability  defined  on  the  basis  of  the  location  of  the  poles  of  the  subsystem  under  consideration. 

In  spite  of  altering  the  location  of  the  threshold  of  stability  of  the  linear  model,  it  was  found  that 
the  non  linear  system  remained  stable  with  gains  much  higher  than  those  found  by  the  root  locus 
method.  The  optimal  gains  obtained  by  using  the  above  procedure  were  of  the  order  of  1G-J0,  much 
smaller  than  gain  values  that  kept  the  non-linear  system  stable.  To  determine  and  confirm  that  the 
location  of  the  threshold  of  stability  resulted  in  very  small  gain  values,  the  root  locus  analysis  was 
repeated  with  the  location  of  the  threshold  of  stability  varied.  A  test  case  with  three  controllers 
spaced  0.5A  apart,  and  with  the  first  controller  0.2<5  downstream  the  stirring  region  was  used  to 
analyse  the  variation  in  gains  with  respect  to  the  variation  in  the  threshold  of  stability.  Figure  33 
indicates  that  as  the  threshold  of  stability  is  moved  farther  away  from  unity,  the  root  locus  paths 
remain  inside  the  defined  stable  region  for  a  longer  time,  and  higher  gains  are  obtained  as  a  result 
of  the  application  of  the  methodology.  This  confirms  the  conclusion  that  the  threshold  of  stability 
of  the  linear  model  does  not  match  well  with  the  region  of  stability  of  the  non-linear  system. 

5.5.3  Accuracy  of  the  linear  model 

To  check  the  accuracy  of  the  linear  model,  the  following  tests  were  carried  out: 

•  The  poles  of  the  linear  model  were  obtained  for  two  different  timesteps  to  see  if  the  model 
is  accurate.  Figure  34  shows  the  location  of  the  poles  of  the  linear  model  in  the  open-loop 
configuration  for  two  cases  in  which  the  system  identification  was  carried  out  at  timesteps 
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Figure  34:  Location  of  the  poles  of  the  linear  model  in  the  open-loop  configuration  for  two  cases 
with  dt  =  0.0005  and  dt  =  0.00025. 

of  dt  =  0.0005  and  dt  —  0.00025.  It  is  seen  from  the  distribution  of  the  poles  of  the 
linear  model  that  there  arc  only  slight  differences  between  the  two  cases,  which  confirms  the 
accuracy  of  the  linear  model. 

•  The  poles  of  the  dynamic  matrix  A  were  obtained  for  a  conventional  DES97  run  to  check  if 
the  distribution  of  the  poles  of  the  linear  mode!  is  similar  to  the  distribution  obtained  with  the 
precursor  simulation  with  stirring  and  with  the  system  in  open-loop  Note  that  in  this  case, 
there  was  no  random  force  being  added  during  the  system  identification  step.  It  was  found 
that  the  poles  of  the  linear  model  so  obtained  had  a  similar  distribution  to  the  earlier  cases. 

•  The  linear  model  was  used  to  predict  future  data  (of  the  resolved  shear  stress)  and  this  data 
obtained  was  compared  with  data  obtained  from  the  non-linear  system.  Figure  35  shows  that 
the  data  obtained  from  the  linear  model  compares  favorably  with  the  data  obtained  from  the 
non-linear  system,  further  confirming  the  accuracy  of  the  linear  model, 

5.5.4  Reduction  of  the  order  of  the  linear  model  using  Proper  Orthogonal  Decomposition 
(POD) 

Proper  Orthogonal  Decomposition  (POD)  is  a  technique  that  can  be  used  to  simplify  a  data  set  to 
a  lower  dimension  for  ease  of  analysis.  This  can  be  accomplished  in  a  manner  that  the  dimension¬ 
ality  of  a  data  set  is  reduced,  but  the  dataset  still  retains  those  characteristics  of  the  data  set  that 
contribute  most  to  its  variance. 

The  main  reasons  for  generating  reduced  models  include: 

•  Lower  order  models  reduce  the  computational  effort  to  obtain  the  desired  results. 

•  Lower  order  models  facilitate  or  ease  controller  design  with  the  resulting  controller  having  a 
simpler  structure  that  is  easier  to  understand  and  parameterize. 
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Figure  35:  Comparison  of  the  output  of  the  linear  model  with  the  output  of  the  system. 

In  the  current  work,  the  previous  section  described  the  problem  with  root  locus  analysis  of 
the  linear  model,  wherein  it  becomes  impossible  to  correctly  define  the  region  of  stability  so  that 
the  threshold  of  stability  of  the  linear  model  matches  with  the  region  of  stability  of  the  non-linear 
system.  One  of  the  ways  to  ameliorate  this  situation  is  to  ensure  that  the  root  locus  analysis  yields 
simpler  paths,  that  can  facilitate  a  more  straightforward  analysis  by  reducing  the  order  of  the  data 
collected  during  the  system  identification  step.  The  application  of  POD  to  reduce  the  order  of  the 
linear  model  is  described  in  detail  below.  The  following  set  of  equations  show'  the  derivation  of  the 
dynamic  matrices  A  and  B  using  the  linear  model. 

X  =  [A  B]  [Xs  rf  (61) 

A'  [Xs  t 1}  =  [A  B][XS  if  \XS  t}]  (62) 

[/l  B]  =  X\XS  7/]  (\XS  r f[Xs  t?])  1  (63) 

From  the  above  equations,  matrices  X  and  X,  contain  the  data  recorded  during  the  system  iden¬ 
tification  step.  The  first  step  is  to  transform  the  given  data  sets  X  and  Xs  into  alternative  data 
sets  of  smaller  dimensions.  The  method  of  Sirovitch  is  used  to  calculate  the  POD  basis  vectors. 
Accordingly,  matrices  B  and  Bx  are  formed  using  X  and  Xs  in  the  following  way: 

B  =  XTX  (64) 

Bs  =  X,TX,  (65) 

Next,  the  eigenvalues  and  eigenvectors  of  the  matrices  B  and  Bs  are  obtained.  The  eigenvector 
corresponding  to  the  highest  eigenvalue  is  used  to  calculate  the  POD  basis  vector.  If  y  and  ys  denote 
the  eigenvectors  corresponding  to  the  highest  eigenvalues,  the  POD  basis  vectors  are  calculated  as, 

V  =  Xy 

^  ,v  A, 
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Figure  36:  Root  locus  plot  using  the  reduced  order  data.  In  this  diagram,  x  represents  the  poles  of 
the  system. 

The  new  data  set  that  is  of  reduced  order  is  recalculated  as, 

Xncw  =  VTX  (68) 

*5l«eu,  =  V,TX'  (69) 

The  new  data  set  of  reduced  order  can  then  be  used  to  calculate  the  dynamic  matrices  A  and  B 
accordingly.  Figure  36  shows  the  root  locus  plot  made  with  the  reduced  order  data.  In  this  case 
the  order  of  the  data  has  been  reduced  to  1,  since  the  eigenvector  that  corresponds  to  the  highest 
eigenvalue  has  been  retained  in  the  data  set.  The  reduction  of  order  can  be  also  accomplished  by 
considering  a  few  of  the  eigenvectors  that  correspond  to  the  largest  eigenvalues.  In  this  case,  a 
single  eigenvector  corresponding  to  the  largest  eigenvalue  was  retained.  To  close  the  loop  a  pole 
is  placed  at  z  =  1.  It  is  now  seen  that  there  are  two  poles  for  the  linear  model  in  figure  36,  one 
pole  at  the  center  of  the  unit  circle  and  the  second  pole  corresponding  to  z  —  1.  With  such  a  pole 
distribution,  the  root  locus  remains  within  the  unit  circle  for  a  longer  time,  resulting  in  higher  gain 
values.  It  was  found  however,  that  the  gain  values  thus  obtained  arc  too  high  and  do  not  ensure  the 
stability  of  the  non-linear  system. 
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5.5.5  Method  based  on  non-linear  system  response 


During  the  next  stages  of  the  investigation,  the  response  of  the  non-linear  system  to  preset  gain 
values  was  measured  during  a  precursor  simulation.  In  the  precursor  simulation,  gain  values  are 
increased  starting  from  zero  until  there  are  numerical  instabilities.  The  gain  values  arc  increased 
on  the  basis  of  a  parameter  d>  that  is  related  to  the  response  of  the  system  where  <p  is  defined  by 


(70), 


M 


Vt.M  ~  Vi-N.M 


Vt-N,M 


N 


Vt 


=  jvDuV)" 


(70) 


The  parameter  <j>(t)  is  measured  at  each  control  plane.  The  controller  gain  is  increased  only  when 
<j>(t)  is  less  than  a  threshold  (5%),  set  to  ensure  numerical  stability.  This  approach  determines 
optimal  gain  values  for  each  control  plane  and  is  similar  to  a  root  locus  method  for  a  non-linear 
system  since  the  gain  values  ensure  stability  of  the  system  based  on  the  response  to  the  actual 
inputs.  To  ensure  the  accuracy  of  the  gains  predicted  using  this  approach,  the  gain  values  must  be 
increased  in  small  increments.  In  the  above  formulation,  the  parameter  M  refers  to  the  number  of 
points  in  the  wall-normal  direction  at  which  the  control  force  is  applied.  Since  the  method  is  based 
upon  the  frequency  of  the  response  measured  via  the  parameter  and  the  gains  are  incremented 
only  when  the  response  has  a  low  frequency,  N  is  a  parameter  similar  to  a  time-averaging  window 
and  was  set  to  a  value  of  100  in  the  simulations  conducted.  With  the  increment  set  to  10-e.  the 
simulation  takes  approximately  15  time  units  to  determine  the  optimal  gain  values.  Despite  the  use 
of  the  small  increment,  the  method  was  not  found  to  be  successful  in  ensuring  the  stability  of  the 
system.  This  necessitated  a  deeper  study  of  the  errors  during  the  control  process  to  possibly  arrive 
at  simpler  alternatives  and  possibly  eliminate  the  necessity  of  a  precursor  simulation. 


5.5.6  Difference  between  PI  controllers  and  integral  controllers 

Here,  results  that  show  that  a  purely  integrator  based  controller  is  as  effective  as  a  PI  controller  in 
increasing  the  resolved  stress  to  target  levels  are  presented.  The  results  shown  below  confirm  this 
conclusion,  by  comparing  the  effect  of  a  controller  using  the  integrator  only  with  other  controllers 
using  different  values  of  Kp. 

Various  values  were  used  for  the  proportional  and  integral  gain  in  a  number  of  simulations  con¬ 
ducted  to  assess  the  effect  of  these  parameters  on  the  system.  Figure  37  shows  plots  that  indicate 
the  effect  of  Kp  and  K;  on  the  resolved  stress  values  at  the  third  control  plane  after  two  time  units. 
Both  plots  indicate  that  the  controller  is  effective  even  if  Kp  =  0  with  a  purely  integral  controller 
acting  on  the  system.  The  figures  also  show  that  the  simulations  with  smaller  values  of  K /  ex¬ 
hibit  the  slowest  response.  These  findings  indicate  that  a  purely  integrator-based  controller  would 
suffice  for  this  system,  obviating  the  need  to  choose  two  constants  (K P  and  A'/).  From  a  consider¬ 
ation  of  steady-state  operation  only,  integrator-based  control  systems  are  preferable  to  proportional 
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Figure  37;  Resolved  stresses  at  the  third  control  plane  for  various  values  of  Kp  and  Kt. 

systems.  However,  it  is  generally  easier  to  achieve  good  transient  behavior  (lower  standard  devia¬ 
tion  of  the  error  about  a  mean)  with  a  proportional-based  control  system  than  an  integrator-based 
control  system.  The  action  of  a  proportional  plus  integral  controller  in  response  to  a  change  in  the 
input  or  external  disturbance  is  initially  similar  to  that  of  a  proportional  controller,  but  as  the  new 
equilibrium  point  is  reached,  the  control  action  becomes  the  same  as  that  of  an  integral  controller. 
A  proportional  plus  integral  controller  combines  the  desirable  transient  characteristics  of  a  propor¬ 
tional  controller  and  the  feature  of  no  steady-state  error.  From  the  point  of  view  of  eliminating  the 
mean  error  in  the  simulation,  an  integral  controller  is  sufficient  in  achieving  the  desired  result. 

5.5.7  Characterization  of  the  terms  contributing  to  the  error 

Figure  38a,b,c  and  d  show  the  effect  of  various  values  of  the  proportional  and  integral  gains  on 
the  time  histories  of  the  mean  of  the  error  52  obtained  at  various  wall-normal  locations.  The 
simulations  were  performed  with  a  fixed  inlet  condition  obtained  from  a  RANS  solution  at  the 
corresponding  Reynolds  number.  The  boundary  layer  stirring  force  ftii  is  then  applied  at  a  single 
stirring  plane  (x/S  —  1.0).  The  simulation  featured  the  use  of  three  control  planes.  Note  that  these 
plots  have  been  obtained  using  data  from  the  third  control  plane  of  the  simulations  in  a  simulation 
conducted  at  a  Reynolds  number  of  5000  and  with  the  single  stirring  plane  spaced  0.2 6  upstream 
of  the  first  control  plane.  The  four  locations  shown  correspond  to  y/S  =  0.0321  (y  +  =  160), 
y/S  =  0.1344  (y+  =  672),  y/S  =  0.3383  (y+  =  1692),  and  y/S  =  0.7177  (y+  =  3589).  The 
figures  show  that  as  time  increases,  the  mean  of  the  errors  goes  towards  zero  at  all  of  the  locations. 
In  the  cases  with  Kp  =  0  (Figures  38b  and  c),  it  can  be  seen  that  the  mean  error  at  locations 
corresponding  to  y/S  =  0.3383  and  y/S  =  0,1344  hardly  decrease  at  the  start  of  the  simulation, 
up  to  a  time  of  approximately  0.5  -  1.0  units,  and  then  decrease  rapidly  tending  towards  zero  as 
time  increases.  In  contrast,  the  cases  with  the  proportional  and  integral  terms  contributing  to  the 
error  show  steady  reduction  in  the  mean  values  of  the  error  over  time  beginning  from  the  start  of 
the  simulation. 
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Figure  39  shows  the  error  histories  at  the  third  control  plane  at  two  wall-normal  locations  cor¬ 
responding  to  y/6  =  0.1344  and  y/S  =  0.3383.  As  shown  in  the  figure,  the  local  errors  are  of 
higher  magnitude  when  a  purely  integral  controller  is  used.  Despite  this,  the  errors  at  these  loca¬ 
tions  have  a  zero  mean  which  eventually  ensure  that  the  resolved  stresses  reach  the  target  values 
and  also  ensure  that  the  quantity  r  does  not  increase  in  magnitude  as  time  increases.  Figure  40 
shows  the  comparison  of  the  rate  of  reduction  of  the  mean  of  the  error  for  all  simulations  at  the 
same  two  locations.  Figure  40a  shows  that  there  are  no  marked  differences  in  rate  of  reduction  of 
the  mean  of  the  error  for  all  the  simulations.  Figure  40b  shows  that  the  addition  of  the  proportional 
term  affects  the  rate  of  reduction  of  the  error  at  the  start  of  the  simulation,  but  as  the  simulation 
proceeds,  the  rate  of  reduction  of  the  errors  is  almost  the  same  for  all  the  cases  shown. 
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Figure  38:  Mean  values  of  the  error  over  time  at  various  wail -normal  locations  at  the  third  control 
plane. 

With  an  integral  controller,  whenever  the  error  has  higher  frequencies,  the  sum  of  the  error  may 
not  take  this  into  account,  especially  when  the  sum  of  the  errors  is  very  large  due  to  a  long  period 
of  low  frequency  errors.  At  the  initiation  of  the  simulation,  the  error  is  more  or  less  constant  and 
shows  low  frequency  content  due  to  which  the  sum  of  the  errors  increases  monoton ically.  When 
the  sum  of  the  errors  becomes  sufficiently  large,  the  errors  begin  to  decrease  and  the  resolved  stress 
moves  towards  the  set  point.  With  larger  values  of  the  integral  gain,  K the  term  7-  becomes  larger 
during  the  initial  stages  of  the  simulation,  resulting  in  a  larger  overshoot  of  the  resolved  stress 
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Figure  39:  Error  histories  at  the  wall-normal  locations  y/6  =  0.3383  and  y/6  =  0.1344  at  the  third 
control  plane. 


Figure  40:  (a)  and  ( b )  Comparison  of  the  rate  of  reduction  of  the  mean  error  over  time  at  wall- 
normal  locations  corresponding  to  y/S  —  0.1344  and  y/6  —  0.3383  at  the  third  control  plane. 
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when  the  error  reaches  zero  for  the  first  time.  In  that  sense,  the  selection  of  an  appropriate  value 
for  the  integral  gain  is  important  as  a  very  large  increase  in  the  error  may  produce  large  variations  in 
the  instantaneous  shear  stress  which  is  undesirable.  However,  the  action  of  the  integral  controller 
ensures  that  the  error  reduces  subsequently  and  moves  in  the  opposite  direction.  This  behavior 
repeats  and  results  in  error  histories  that  are  characterized  by  large  modulations  in  the  error  for  a 
very  short  period  of  time,  followed  by  long  periods  of  low  frequency  errors.  Despite  this  behavior, 
the  average  error  goes  to  zero  after  some  time  and  undesirable  transient  effects  (primarily  on  the 
wall  shear  stress)  can  be  avoided  by  the  use  of  smaller  values  of  integral  gains  (typically  values 
ranging  from  0.2  to  0,5  yielded  the  desired  results).  From  this  perspective,  an  integral  controller 
is  sufficient  to  provide  the  desired  results  although  the  gain  values  need  to  be  selected  carefully. 
With  the  use  of  a  PI  controller,  the  r  term  attains  targe  values  at  the  beginning  of  a  simulation,  The 
higher  frequency  components  of  the  error  are  taken  into  account  in  determining  the  control  force 
due  to  the  presence  of  the  proportional  term.  The  integral  term  is  still  present  to  account  for  the  low 
frequency  components  of  the  error.  The  value  of  Kp  plays  a  role  in  reducing  the  high  amplitude 
modulations  of  the  error  that  occurs  when  an  integral  controller  is  used.  The  main  differences  in 
the  error  histories  of  the  PI  controller  and  integral  controller  are  contrasted  in  Figure  41  which 
considers  the  error  histories  at  all  three  control  planes  for  two  simulations  with  the  same  integral 
gain  of  Kj  =  0.5  and  Kp  =  0.02.  Note  that  these  simulations  were  carried  out  at  a  Reynolds 
number  Rer  =  5000  and  with  a  single  stirring  plane  and  three  controllers  spaced  0.55  apart.  The 
errors  at  wall-normal  locations  corresponding  to  y/5  =  0  1344  (y+  =  672)  and  y/5  —  0,3383 
(y+  =  1692)  are  shown  in  the  figure.  The  response  of  the  system  shows  smaller  overshoots 
when  PI  controllers  are  used.  The  use  of  integral  controllers  requires  the  selection  of  just  a  single 
parameter  -  the  integral  gain  and  is  therefore  advantageous. 

Further  investigations  that  were  conducted  to  understand  the  behavior  of  the  error  at  the  control 
planes  included: 

•  Effect  of  the  stirring  profile  on  the  error  histories.  With  a  modified  stirring  profile  (c.f..  Fig¬ 
ure  55)  the  histories  of  the  error  and  the  quantity  r  showed  no  changes  in  terms  of  behavior 
(in  terms  of  the  history  of  the  errors  overtime).  Note  that  for  the  comparison,  the  simulations 
were  conducted  with  values  of  Kp  and  K\  corresponding  to  0.2  and  0.5  respectively. 

■  Effect  of  the  spacing  between  the  end  of  the  stirring  region  and  the  first  control  plane  on 
the  error  histories.  With  the  increase  in  the  streamwise  distance  between  the  end  of  the 
stirring  region  and  the  first  control  plane  to  5,  the  error  histories  at  the  control  planes  did  now 
show  any  marked  differences  when  compared  with  the  cases  with  0.25  streamwise  distance 
between  the  end  of  the  stirring  region  and  first  control  plane. 
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Figure  41 :  Comparison  of  error  histories  at  two  wall-normal  locations  at  all  three  control  planes 
when  PI  controllers  and  integral  controllers  are  used. 
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5.5.8  Non-linear  PI  controller 


The  amplitudes  of  the  error  modulations  can  also  be  reduced  by  using  a  non-linear  PI  controller 
with  the  proportional  gain  KP  given  by  the  rate  of  change  of  error  (de/dt).  This  is  based  on 
the  premise  that  when  the  error  has  a  low  frequency  (and  is  not  changing  much),  the  integral 
term  is  sufficient  to  provide  the  desired  control  as  de/dt  =  0.  When  the  error  exhibits  a  high 
frequency,  de/dt  can  be  large,  depending  upon  the  frequency,  resulting  in  greater  proportional 
action  that  dominates  over  the  integral  action.  This  non-linear  controller  can  be  interpreted  as 
a  hybrid  between  a  PI  controller  and  an  Integral-Differential  (ID)  controller,  as  the  formula  that 
expresses  the  control  action  can  be  written  as: 


(71) 


When  de/dt  =constant,  the  controller  behaves  as  a  PI  controller  and  when  e  is  nearly  constant,  the 
controller  becomes  an  approximation  of  a  ID  controller.  When  e  is  exactly  equal  to  a  constant,  the 
controller  behaves  as  an  integral  controller. 

In  this  way,  one  could  obtain  a  significant  (i.e.,  different  from  zero)  control  action  when  the 
error  is  nonzero  and  when  there  are  high  frequency  oscillations  of  the  error.  Furthermore,  if  one 
examines  the  term  e(de/dt),  it  is  apparent  that  the  amplitude  of  the  control  action  always  varies 
in  a  way  that  prevents  the  controller  gain  from  becoming  infinite.  Figure  42  shows  the  reduction 
in  the  amplitude  of  the  error  modulations  at  the  wall-normal  location,  y/S  —  0.3383,  due  to  the 
incorporation  of  the  proportional  gain  KP  =  de/dt.  Subsequent  simulations  were  carried  out 
to  determine  the  sensitivity  of  the  choice  of  KP  on  the  errors.  In  these  simulations,  values  of 
KP  corresponding  to  2de/dt,  Ade/dl  and  8 de/dt  were  used.  Figure  43  shows  the  comparison  of 
the  error  histories  at  a  wall-normal  location  corresponding  to  y/8  —  0.3383  for  two  cases  with 
KP  =  2 de/dt  and  KP  -  de/dt.  With  KP  =  2 de/dt,  the  amplitude  of  the  error  modulations 
are  slightly  smaller.  Figure  44  shows  a  similar  comparison  for  two  cases  with  K P  —  4 de/dt  and 
KP  —  Sde/dt.  With  KP  —  8  de/dt,  the  amplitude  of  the  error  modulations  are  the  smallest.  These 
cases  serve  to  illustrate  that  the  choice  of  the  factor  that  multiplies  de/dt  is  not  critical  especially 
in  the  absence  of  any  special  advantage  in  reducing  the  amplitude  of  the  error  modulations  further 

5.5.9  Effect  of  time-averaging  in  the  control  formulation 

Equation  (52)  is  used  to  calculate  the  error  at  each  control  plane.  The  equation  show's  that  the  error 
is  the  difference  between  the  set-point  and  the  average  value  of  resolved  shear  stress.  The  type 
of  averaging  employed  is  a  moving  time-average  and  simulations  were  conducted  to  assess  the 
cfTect  of  the  width  of  the  time-averaging  envelope  on  the  results.  Simulations  were  conducted  with 
Ki  =  0,5  and  the  width  of  the  time-averaging  envelope  varied  from  1  timcstcp(no  time-averaging) 
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A«  yrb  =  0  3303.  Kp  =  0  and  0  5 


Figure  42:  Reduction  in  the  amplitude  of  the  error  modulations  due  to  the  action  of  the  non-linear 
PI  controller. 
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Figure  43:  Comparison  of  the  error  histories  for  two  cases  corresponding  to  Kp  =  2de/dt  and 
K /  =  0.5;  Kp  =  de/dt  and  K{  =  0.5. 
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At  y/&  =  0  3383,  kp  =  4  (tefdt,  R,  =  C  5 


Figure  44:  Comparison  of  the  error  histories  for  two  cases  corresponding  to  Kp  =  Ade/dt  and 
K,  =  0.5;  KP  =  Sde/dt  and  K,  =  0.5. 

to  200  timestcps  (corresponding  to  0.05  time  units).  Plots  of  the  Power  Spectral  Density  (PSD)  of 
the  error  signals  were  analysed  to  determine  the  effect  of  time-averaging  on  the  frequency  content 
seen  by  the  Integrators. 

Figure  45a  shows  the  PSD  of  the  error  signal  at  a  location  of  y/d  —  0.3383  at  the  third  control 
plane.  The  figure  shows  that  an  increase  in  the  width  of  the  time-averaging  envelope  results  in  a 
smaller  peak  frequency  of  the  error  seen  by  the  Integrator.  For  example,  with  a  time-averaging 
envelope  of  0.05  time  units,  the  highest  frequency  errors  seen  by  the  Integrator  correspond  to 
timescales  in  the  flow  of  about  0.0333  time-units  (approximately  133  limesteps).  This  implies  that 
all  fluctuations  that  correspond  to  smaller  timescales  are  not  seen  by  the  Integrator  and  are  atten¬ 
uated  as  a  result  of  the  averaging  process.  However,  with  no  time-averaging,  the  peak  frequency 
of  the  error  seen  by  the  Integrator  corresponds  to  a  flow  timescale  of  approximately  0.0077  time 
units  (about  30  timesteps).  Higher  frequencies  seen  by  the  Integrator  imply  better  control  of  the 
resolved  shear  stress.  With  the  current  value  oflntegrator  gain  of  0.5,  increasing  the  width  of  the 
lime-averaging  envelope  to  0.10  and  0.15  time  units  resulted  in  numerical  instabilities.  Simulations 
were  also  conducted  with  Kj  =  0.1  and  with  larger  time-averaging  envelopes  of  1000,  2000  and 
3000  steps  (corresponding  to  0.25,  0.5  and  0-75  lime  units).  Figure  45b  shows  a  PSD  plot  of  the 
errors  at  a  wall-normal  location  of  y/d  —  0.3383  at  the  third  control  plane.  The  peak  frequency 
of  the  error  seen  by  the  Integrator  corresponds  to  timescales  in  the  flow  of  about  0.0667  time  units 
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(about  265  timesteps).  This  implies  that  for  errors  over  all  scales  of  the  flow  that  have  a  timescale 
lower  than  0.0667  time-units  arc  not  taken  into  consideration  by  the  Integrator  and  arc  attenuated. 


Figure  45:  Power  spectra  of  the  errors  with  different  averaging  envelopes. 


5.5.10  Adaptive  control  methodology  to  specify  the  controller  gains 

The  previous  simulations  described  until  now  showed  the  absence  of  major  differences  whether 
a  PI  controller  or  an  integral  controller  is  used,  provided  the  integral  gain  is  selected  carefully. 
To  further  improve  the  selection  methodology,  integral  controllers  were  used  with  the  gain  being 
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varied  adaptively  during  the  simulation.  Initial  investigations  of  the  adaptive  control  methodology 
were  carried  out  with  the  gain  being  varied  at  each  wall-normal  location  based  on  the  error  at  the 
wall-normal  location.  The  formulation  used  for  the  control  input  is. 


r  —  Kj  I  edt 

(72) 

Jo 

l<i,, +,  =  Ki,,  ~  edt 

(73) 

i,i+i  =  Ki,,  +  edt 

(74) 

The  sign  of  the  error  plays  an  important  role  in  determining  the  gain  values  in  the  above  equa¬ 
tions.  At  the  initiation  of  the  simulation,  for  the  bottom  half  of  the  channel,  the  error  is  negative 
(target  resolved  shear  stress  value  is  negative).  Equation  (73)  is  used  in  this  case  while  (74)  is  used 
for  the  upper  half  of  the  channel.  At  the  start  of  the  simulation,  the  error  is  large  due  to  which 
the  gain  increases  by  larger  amounts.  As  the  error  starts  decreasing  due  to  the  action  of  the  inte¬ 
gral  controller,  the  increase  in  gains  become  smaller.  When  the  error  reaches  zero  and  there  is  an 
overshoot,  the  gain  begins  to  reduce.  The  initial  value  of  the  gain  is  the  only  choice  made  prior 
to  the  simulation.  Figure  46  shows  the  error  histories  at  wall-normal  locations  corresponding  to 
y/5  —  0.1344  and  y/6  =  0.3383  using  the  adaptive  control  methodology  with  an  initial  value  of 
K,  =  0.2. 


Al*J6=  0  1344  K*0?  *  n  dl 


Figure  46:  Error  histories  at  two  wall-normal  locations  using  the  adaptive  control  methodology 
with  Kj  —  0.2 
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A  more  standard  form  of  an  adaptive  controller  can  be  used  which  is  described  below.  The  gain 
is  varied  at  each  wall-normal  location  based  on  the  error  at  the  corresponding  location.  The  method 
is  based  on  the  MIT  rule  which  is  an  adaptive  control  methodology  proposed  by  Whitaker  (75). 
The  formulation  used  for  the  control  input  consists  of  a  cost  function,  defined  below  in  (75).  The 
integral  gain  is  calculated  using  (76),  with  the  system  in  open-loop  at  the  start  of  the  simulation 

(JO  =  0). 

In  (77),  in  order  to  prevent  the  gains  from  becoming  too  large,  a  modification  is  introduced 
with  the  gain  being  calculated  from  (78)  every  time  the  value  of  the  term  de/dKin_t  \  is  larger 
than  unity.  It  must  be  noted  that  the  sign  preceding  tj  is  negative  only  for  the  lower  half  of  the 
channel.  For  the  upper  half  of  the  channel,  the  sign  preceding  tj  is  positive.  These  modifications 
have  been  made  to  ensure  that  the  gain  of  the  system  increases  in  the  positive  direction  from  the 
start  of  the  simulation.  During  the  simulation,  by  virtue  of  the  values  of  the  resolved  shear  stress 
and  errors  at  the  control  planes,  if  the  gains  become  negative,  the  gain  values  arc  set  to  zero  so  that 
the  system  gain  values  are  always  positive.  In  this  method,  the  only  choice  to  be  made  by  the  user 
in  such  a  simulation  is  the  value  of  tj  that  ensures  numerical  stability. 


C  = 


(75) 


A'/„  =  JO..,  -  V 

Kin  =  KIn_ ,  -  T/tn 


dC 


dK 

den 


dK,„_ 


(76) 

(77) 


K,  „  =  K,n_%  -  tj  dC  (78) 

A  few  values  of  tj  were  tested  and  the  results  are  shown  in  Figure  47.  Figure  47  shows  the  mean 
of  the  error  at  the  last  control  plane  at  wall-normal  locations  corresponding  to  y/6  =  0.1344  and 
y/6  =  0.3383  using  the  adaptive  control  methodology  for  values  of  tj  =  1,  3,  and  10.  The  figure 
shows  that  the  mean  of  the  errors  goes  to  zero  faster  starting  from  the  system  being  in  the  open 
loop  when  tj  —  10.  With  the  lower  values  of  tj,  it  can  be  seen  that  the  mean  of  the  errors  moves 
very  slowly  towards  zero.  At  larger  values  of  tj  greater  than  10,  the  change  in  the  gain  values  at 
successive  timesteps  becomes  too  large  resulting  in  the  generation  of  numerical  instabilities. 


5.6  Flow  visualizations 

Figure  48  shows  a  comparison  of  the  streamwise  velocity  contours  at  three  wall-normal  locations 
(corresponding  to  y+  =161,  672  and  1692)  for  the  fine  grid  case  (consisting  of  129  x  75  x  65 
points  in  the  streamwise,  wall-normal  and  spanwise  directions  respectively).  Figures  49  shows 
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Figure  47:  (a)  Mean  of  the  error  at  two  wall-normal  locations  at  the  third  control  plane  using  the 
adaptive  control  methodology  with  jj  —  1;  (b)  Mean  of  the  error  at  a  wall-normal  y/5  —  0,1344  at 
the  third  control  plane  using  the  adaptive  control  methodology  with  //  =  1  (Zoomed  in  view);  (c) 
Mean  of  the  error  at  two  wall-normal  locations  at  the  third  control  plane  using  the  adaptive  control 
methodology  with  rj  —  3;  (d)  Mean  of  the  error  at  two  wall-normal  locations  at  the  third  control 
plane  using  the  adaptive  control  methodology  with  tj  =  10. 
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the  contours  of  the  spanwise  velocity  at  the  same  three  wall-normal  locations.  Figure  50  shows 
comparisons  of  the  streamwise  vorticity  at  three  wall-normal  locations  (corresponding  to  i/4  =27, 
161  and  1692).  All  contours  show  evidence  of  turbulent  structures  downstream  the  control  planes. 
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Figure  48:  Streamwise  velocity  contours  at  wall-normal  locations  corresponding  to  y/6  —  0.0321 
(y4  =  161),  y/6  =  0.1344  (y+=672)  and  yfS  =  0.3383  (y+  =  1692). 


5,7  Resolved  stress  and  wall  stress 

Figure  51a  shows  the  resolved  shear  stress  at  the  last  control  plane  from  a  simulation  performed 
on  a  channel  with  a  streamwise  extent  of  4tt,  with  a  single  stirring  plane  at  z/6  —  1.0  and  five 
controllers  spaced  6  apart.  The  simulation  was  performed  on  the  fine  grid  consisting  of  129  x 
75  x  65  points  in  the  streamwise,  wall-normal  and  spanwise  directions  respectively.  The  plot 
shows  that  the  peak  resolved  shear  stress  reaches  within  10%  of  the  set  point  by  about  two  time 
units.  Figure  51b  shows  the  resolved  shear  stress  measured  at  streamwise  locations  3c>  and  5 6 
downstream  the  last  control  plane.  It  can  be  seen  that  the  peak  values  of  the  resolved  shear  stress 
at  these  locations  are  very  close  to  the  peak  resolved  stress  of  the  set  point  (target  profile),  even 
though  the  profiles  seem  to  be  shifted  to  the  left  with  the  peak  resolved  stress  at  the  locations 
downstream  occurring  at  a  location  corresponding  to  y/6  =  0.18  as  compared  to  a  location  of 
y/6  =  0.3  in  the  target  profile. 
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Figure  49:  Spanwise  velocity  contours  at  wall-normal  locations  corresponding  to  y/S  =  0.0321 
(y+  =  161),  y/S  =  0.1344  (y+=672)  and  y/S  =  0.3383  (y+  =  1692). 
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Figure  50:  Streamwise  vorticity  contours  at  wall-normal  locations  corresponding  to  i/+=27,16t 
and  1692. 
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Figure  51c  shows  the  resolved  shear  stress  measured  at  streamwise  locations  3(5  and  5<5  down¬ 
stream  of  the  last  control  plane  from  a  simulation  conducted  on  the  coarse  grid  consisting  of 
129  x  75  x  33  points  in  the  streamwise,  wall-normal,  and  spanwise  directions,  respectively.  As 
in  the  previous  case,  this  simulation  was  performed  on  a  channel  with  a  streamwise  extent  of  At, 
and  with  a  single  stirring  plane  at  x/t 5  —  1.0  and  five  controllers  <5  apart.  The  plot  shows  that  the 
peak  values  of  the  resolved  shear  stress  at  the  locations  3(5  and  5(5  downstream  of  the  last  control 
plane  are  substantially  less  than  the  peak  values  of  the  resolved  shear  stress  of  the  target  profile  by 
approximately  15%  and  25%,  respectively.  This  reduction  in  the  peak  values  of  the  resolved  shear 
stress  downstream  of  the  control  plane  imply  an  under-prediction  in  the  wall  stress. 

Further  investigations  revealed  that  the  under-prediction  in  the  wall  stress  is  a  result  of  the 
use  of  grids  with  insufficient  density  and  the  number  of  control  planes  used  for  the  simulation. 
Figure  52  shows  a  comparison  of  the  wall  shear  stress  for  the  coarse  grid  and  fine  grid  cases 
(coarse  grid  consisting  of  129  x  75  x  33  points  and  fine  grid  consisting  of  129  x  75  x  65  points 
in  the  streamwise,  wall-normal,  and  spanwise  directions  respectively).  The  channel  was  long 
in  the  streamwise  direction  and  the  stirring  force  was  applied  at  a  single  plane,  x/5  —  1.0  for 
both  simulations.  The  finer  grid  had  65  points  in  the  spanwise  direction,  resulting  in  =  0.0491 
(roughly  20  points  per  boundary  layer  thickness).  In  this  way,  the  location  of  the  RANS-LES 
interface  for  the  finer  grid  is  located  at  a  distance  of  y/S  =  0.0318,  corresponding  to  a  of 
approximately  160.  The  simulations  were  conducted  at  a  Reynolds  number  of  5000.  The  figure 
shows  that  in  the  case  of  the  fine  grid,  there  is  a  good  recovery  of  the  wall  shear  stress  to  within 
7  —  8%  of  the  correct  levels,  with  relatively  small  decreases  in  wall -shear  stress  downstream  of  the 
control  planes.  The  coarse  grid  levels  are  low  by  about  20%  close  to  the  exit  of  the  channel  due  to 
a  steady  decay  in  the  wall  stress  as  the  flow  moves  downstream. 

Figure  53  shows  comparisons  of  the  wall  shear  stress  for  three  fine-grid  cases  with  variations 
in  the  number  of  control  planes  and  the  spacing  between  them.  All  the  runs  were  performed  on  a 
domain  with  a  streamwise  extent  of  At  and  a  grid  with  20  points  per  boundary  layer  thickness  in 
the  spanwise  direction.  It  can  be  seen  that  with  5  controllers  spaced  <5/2,  the  wall  shear  stress  drops 
steadily  downstream  of  the  control  planes  with  an  error  of  approximately  11%  at  the  exit  of  the 
domain.  Despite  the  drop  in  the  wall  shear  stress  being  less  marked  as  the  coarse  grid  result,  it  is  an 
undesirable  effect.  With  5  controllers  spaced  6  or  with  10  controllers  spaced  c5/2,  the  reduction  in 
the  wall  shear  stress  is  lower  with  the  error  close  to  8%  at  the  exit  of  the  domain.  Figure  54  shows 
distributions  of  the  resolved  shear  stress  at  locations  3<5  and  5(5  downstream  of  the  last  control 
plane.  For  the  case  with  5  controllers  spaced  <5/2  apart,  the  peak  values  of  the  resolved  shear  stress 
are  substantially  lower  than  the  peak  value  of  the  target  profile  indicating  that  there  is  a  decay  of 
turbulent  content  downstream  of  the  control  planes.  The  resolved  shear  stress  for  the  other  two 
cases  (with  10  controllers  spaced  5/2  and  5  controllers  spaced  <5)  do  not  show  a  substantial  drop  in 
the  peak  value  as  compared  to  the  target  profile  and  this  translates  to  a  better  prediction  of  the  wall 
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Figure  5 1 :  Resolved  shear  stress  (a)  with  the  fine  grid  ( 1 29  x  75  x  65  points  in  the  streamwise.  wall- 
nonnal  and  spanwise  directions)  at  the  last  control  plane;  (5)  with  the  fine  grid  (129  x  75  x  65  points 
in  the  streamwise,  wall-normal  and  spanwise  directions)  at  locations  corresponding  to  35  and  55" 
downstream  the  last  control  plane;  (c)  with  the  coarse  grid  (129  x  75  x  33  points  in  the  streamwise, 
wall-normal  and  spanwise  directions)  at  locations  corresponding  to  35  and  55  downstream  the  last 
control  plane. 


5.8  Influence  of  the  stirring  force 

The  influence  of  the  stirring  force  on  the  control  methodology  was  investigated  using  simulations 
performed  at  Rer  —  400.  Computations  were  performed  using  different  profiles  for  referred 
to  as  “scaling  profile  1”  and  “scaling  profile  2”  in  Figure  55.  Scaling  profile  2  was  applied  to 
distribute  the  force  in  the  streamwise  direction  and  the  same  profile  with  a  magnitude  half  that 
used  in  the  streamwise  direction  was  employed  to  scale  the  wall-normal  and  spanwise  components 
of  the  force.  This  anisotropic  treatment  was  used  since  the  streamwise  fluctuations  are  larger  than 
the  wall-normal  and  spanwise  components  in  boundary  layers.  The  results  obtained  using  scaling 
profile  2  are  compared  with  those  obtained  using  scaling  profile  1  (with  the  same  profile  being 
applied  in  all  three  directions). 

Figure  56  and  Figure  57  show  the  nns  levels  of  the  velocity  fluctuations  at  locations  15  and  25 
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Figure  52:  Comparison  of  the  wall  shear  stress  distribution  for  the  coarse  grid  and  fine  grid  cases. 


Figure  53:  Comparison  of  the  wall  shear  stress  distribution  for  three  cases  with  different  number 
and  placement  of  controllers  using  the  fine  grid. 
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Figure  54:  Resolved  shear  stress  (a)  with  the  fine  grid  (129  x  75  x  65  points  in  the  streannvise, 
wall-normal  and  spanwise  directions)  at  a  location  3£  downstream  the  last  control  plane;  (6)  with 
the  fine  grid  (129  x  75  x  65  points  in  the  streamwise,  wall-normal  and  spanwise  directions)  at  a 
location  5«5  downstream  the  last  control  plane. 

downstream  of  the  last  controller  using  the  two  profiles.  The  figure  shows  that  the  velocity  fluctu¬ 
ations  in  this  case  with  stirring  profile  2  exhibit  better  agreement  with  the  rms  velocity  fluctuations 
from  a  fully-developed  channel  flow  simulation  at  the  same  Reynolds  number.  This  indicates  that 
the  form  of  the  stirring  force  plays  a  role  in  ensuring  that  the  flow  downstream  of  the  controller 
remains  realistic. 


Figure  55:  Scaling  profiles  of  the  stirring  force. 


5.9  Synthetic  turbulence  generation  with  controlled  forcing 

The  synthetic  turbulence  generation  method  of  Batten  el  al.  (69)  was  investigated  and  used  to 
provide  a  three-dimensional,  unsteady  velocity  field  at  the  inflow  plane.  The  synthetic  turbulence 
method  was  used  at  the  inlet  (instead  of  boundary  layer  stirring  to  generate  three-dimensional 
content)  along  with  controllers  downstream,  that  amplify  the  fluctuations  to  obtain  target  resolved 
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Figure  56:  RMS  velocities  obtained  using  scaling  profile  1  for  the  stirring  force,  (o)  profiles  Id' 
downstream  of  the  last  controller;  ( b )  profiles  2 5  downstream  of  the  last  controller. 


Figure  57:  RMS  velocities  obtained  using  scaling  profile  2  for  the  stirring  force,  (a)  RMS  velocities 
at  a  location  corresponding  to  1<5  downstream  of  the  last  controller;  (6)  RMS  velocities  at  a  location 
corresponding  to  2 5  downstream  of  the  last  controller, 

stress  levels  at  the  control  planes.  The  results  showed  that  synthetic  turbulence  was  successful  in 
generating  fluctuations  that  are  more  realistic  than  boundary  layer  stirring,  as  a  result  of  which  the 
root  mean  square  (nns)  values  of  the  velocities  showed  reasonable  agreement  with  profiles  from  a 
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periodic  simulation  using  DES97.  The  controllers  still  have  a  role  to  play,  and  help  to  ensure  that 
the  levels  of  resolved  stresses  reach  the  the  target  levels  quickly.  Figure  58  shows  the  rms  values 
of  the  streamwise  velocity  at  various  locations  downstream  of  the  inlet.  Figure  59  shows  the  rms 
values  of  the  wall-normal  velocity  at  various  streamwise  locations  downstream  of  the  inlet.  It  can 
be  seen  that  the  profiles  adjust  themselves  within  7-8(5  downstream  of  the  inlet. 


Figure  58:  RMS  values  of  the  streamwise  velocity  at  various  streamwise  locations  (Synthetic 
turbulence  with  controlled  forcing). 
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Figure  59:  RMS  values  of  the  wall-normal  velocity  at  various  streamwise  locations  (Synthetic 
turbulence  with  controlled  forcing). 
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6  Summary 

Prediction  of  complex  flows  that  include  adverse  pressure  gradients,  streamline  curvature  and 
boundary  layer  separation  remain  among  the  most  challenging  for  turbulence  simulation  strate¬ 
gics.  Among  the  various  hierarchies  of  simulation  strategies,  RANS  approaches  remain  the  most 
widely  applied  simulation  techinque  for  complex  flows.  However,  this  approach  often  produces 
unsatisfactory  results  for  flows  far-removed  from  the  thin  shear  layers  used  to  calibrate  the  un¬ 
derlying  RANS  models.  For  example,  significant  effects  of  separation  in  flows  that  are  massively 
separated  are  not  easily  predicted  using  RANS  approaches.  The  strong  reliance  on  empircism 
has  motivated  intense  development  of  other  approaches  that  provide  more  robustness  and  fidelity. 
Strategics  such  as  LES  are  attractive  from  the  point  of  view  that  they  minimize  the  empirical  input 
in  the  models  and  are  capable  of  generating  solutions  with  rich  turbulence  physics.  Despite  this 
benefit,  LES  carries  a  prohibitive  computational  cost  for  resolving  boundary  layer  turbulence  at 
high  Reynolds  numbers  with  grid  point  estimates  that  scale  as  Re1 *  (25).  This  provided  a  strong 
incentive  to  merge  the  two  techniques  in  such  a  way  that  the  advantages  of  both  strategics  are 
preserved.  Among  the  most  popular  of  all  the  hybrid  RANS/LES  techniques  is  DES  (26)  which 
has  been  actively  applied  to  simulate  complex  flowfields  with  massive  separation.  Even  though  the 
method  was  initially  proposed  by  Spalart  et  at.  (26)  as  a  numerically  feasible  strategy  to  predict 
massively  separated  flows,  interest  has  been  generated  to  extend  the  method  to  other  classes  ot 
flows. 

The  use  of  DES  to  predict  attached  flows  or  flows  with  shallow  separation  is  more  challenging 
due  to  the  difficulty  in  the  generation  and  sustenance  of  instabilities.  The  over-arching  goal  ol  the 
investigations  was  to  understand  this  aspect  of  DES  and  understand  the  method  better  so  that  it  can 
be  extended  to  all  classes  of  flows.  The  initial  stages  of  the  project  were  focused  on  improving  the 
overall  understanding  of  DES  for  flows  with  thick  boundary  layers  and/or  shallow  separation.  It 
was  found  that  finer  grids  in  the  streamwise  and  spanwise  directions  caused  the  RANS-LES  inter¬ 
face  to  lie  well  within  the  boundary  layer,  teading  to  a  reduction  of  modeled  stresses  in  the  RANS 
region  of  the  simulation  with  resolved  content  lacking  in  the  LES  region  due  to  insufficient  grid 
densities  in  the  LES  region.  This  further  led  to  a  reduction  in  skin  friction  resulting  in  premature 
separation  of  the  flow.  The  problems  that  occur  in  these  regimes  due  to  the  use  of  ambiguous  grids 
were  eliminated  by  the  use  of  DDES,  the  newer  version  of  the  model.  The  new  version  was  used 
to  simulate  attached  flow  over  a  flat  plate  at  a  momentum  thickness  based  Reynolds  number  of 
9000.  The  new  version  was  then  tested  on  the  flow  over  an  Aerospatiale  A -airfoil  at  an  angle  ol 
attack  of  13.3°  and  at  a  chord-based  Reynolds  number  of  Re  =  2.1  x  106.  To  ensure  that  the  new 
version  does  not  degrade  predictions  in  massively  separated  flows,  (he  flow  over  a  circular  cylinder 
at  Reynolds  numbers  corresponding  to  1.4  x  10s  and  8  x  10°  was  also  analyzed.  These  simulations 
reinforced  the  effectiveness  of  DDES  in  ensuring  that  RANS  was  maintained  within  the  boundary 
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layer,  leading  to  correct  levels  of  modeled  stresses  in  the  hybrid  simulation  as  compared  to  a  RANS 
simulation. 

Another  important  issue  in  the  simulation  of  realistic  flows  using  DES  involves  the  generation 
of  instabilities  in  a  numerical  simulation  in  flowfields  that  do  not  include  any  slope  discontinu¬ 
ity  (presence  of  slope  discontinuities  would  lead  to  the  quick  growth  of  any  generated  instabili¬ 
ties),  This  becomes  an  important  issue  in  simulations  of  spatially  developing  flows  in  a  channel 
or  a  backward  facing  step,  where  the  generation  of  inflow  conditions  that  enable  the  growth  and 
sustenance  of  turbulence  is  paramount.  In  the  past,  researchers  have  developed  various  methods 
of  specifying  inflow  boundary  conditions  meant  to  be  used  in  conjunction  with  wall-layer  models 
for  these  flowfields.  The  current  work  provides  a  method  to  generate  instabilities  in  the  transition 
from  a  RANS  treatment  to  an  LES  treatment  along  the  main  flow  direction  in  a  DES,  which  can 
result  in  applications  due  to  stream  wise  grid  refinement  as  necessitated  by  a  change  in  geometry  or 
by  design  into  the  grid  the  finer  spacing  needed  to  resolve  fluctuations  in  the  hopes  of  improving 
accuracy. 

The  main  focus  of  the  project  involved  the  application  of  a  controls-based  methodology  to 
develop  velocity  fluctuations  for  boundary  layer  simulations  and  testing  of  the  methodology  using 
DES  of  turbulent  channel  flow.  The  method  employs  a  series  of  control  planes  in  which  a  body 
force  is  applied  to  amplify  velocity  fluctuations  seeded  upstream.  The  control  force  is  added  to  the 
wall-normal  momentum  equations  at  the  control  planes,  based  on  a  prescribed  target  profile  tor  the 
resolved  shear  stress.  The  main  focus  of  the  present  manuscript  was  the  design  of  the  controllers 
and  simplifying  the  procedure.  The  design  of  controllers  using  a  linear  model  was  unsuccessful 
due  to  the  inherent  non-linearities  of  the  system  that  ensured  that  the  threshold  of  stability  of 
the  linear  model  did  not  match  well  with  the  region  of  stability  of  the  non-linear  system.  The 
control  methodology  based  on  the  measured  response  of  the  non-linear  system  was  also  found  to 
be  unsuccessful  in  predicting  gain  values  that  ensure  numerical  stability. 

A  deeper  analysis  of  the  behavior  of  the  errors  with  the  use  of  PI  and  integral  controllers  re¬ 
vealed  that  there  are  no  specific  performance  advantages  due  to  the  use  of  PI  controllers  over 
integral  controllers.  The  use  of  PI  controllers  involves  the  specification  of  two  parameters  -  the 
Proportional  gain  (KP)  and  the  integral  gain  (Kr).  Simulations  conducted  revealed  that  the  spec¬ 
ification  of  gains  that  ensure  numerical  stability  was  very  subtle  and  complicated.  Using  integral 
controllers  with  the  gain  of  the  controllers  increased  adaptively  is  a  simple  alternative  that  ensures 
that  the  desired  results  are  obtained.  The  adaptive  control  scheme  developed  and  tested  worked 
satisfactorily  and  yielded  the  desired  results  without  the  use  of  saturations  to  prevent  the  forces 
from  being  applied  based  on  predefined  criteria.  The  method  was  tound  to  work  satisfactorily  with 
the  resolved  stress  reaching  target  levels  at  the  control  planes,  with  the  peak  value  of  the  resolved 
stress  being  within  10%  of  the  peak  value  of  the  target  profile  at  locations  downstream  of  the  con¬ 
trol  plane  provided  sufficient  number  of  controllers  and  finer  grids  were  used.  With  coarse  grids, 
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there  seemed  be  a  larger  decay  in  the  peak  values  of  the  resolved  stress  at  locations  even  2  S  down¬ 
stream  the  control  plane.  This  also  resulted  in  a  large  error  of  about  20%  in  the  wall  stress  at  the 
end  of  the  domain  (at  a  streamwise  distance  of  8 S  from  the  Iasi  control  plane).  Subsequent  simula¬ 
tions  revealed  that  the  decay  of  the  resolved  stress  and  wall  stress  downstream  the  control  planes 
is  largely  dependent  upon  the  mesh  resolution  used  while  applying  the  method.  It  was  found  that 
with  20  points  per  boundary  layer  thickness  in  the  spanwise  direction,  the  wall  shear  stress  levels 
are  within  8%  of  the  RAMS  levels  for  a  simulation  with  a  domain  streamwise  extent  of  4tt  and 
controllers  occupying  30%  of  streamwise  extent  of  the  domain. 
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Appendix 

Cohalt 

Coball  is  a  parallel,  implicit,  unstructured  flow  solver  developed  by  the  Computational  Sciences 
Branch  of  the  Air  Force  Research  Laboratory.  The  algorithm  ot  Cobalt  is  based  on  Godunov’s 
first-order  accurate,  exact  Riemann  method  (59).  The  solver  is  second  order  accurate  temporally 
and  spatially,  and  the  time  stepping  is  implicit.  The  solver  adopts  a  cell-centered,  finite-volume 
approach.  Cobalt  is  compatible  with  grids  of  arbitrary  cell  shapes  in  two  or  three  dimensions. 
Cobalt  was  chosen  over  other  commercial  codes,  as  it  is  compatible  with  both  DES  and  URANS 
models.  The  computation  of  the  flow  solution  consists  of  five  fundamental  tasks: 

•  Construction  of  initial  conditions  for  the  Riemann  problem  at  any  given  face 

•  Solution  of  this  Riemann  problem 

•  Construction  of  viscous  fluxes  at  any  given  face 


•  Time  integration 


•  Boundary  conditions 


The  Euler  and  Navier-Stokes  equations  are  solved  in  an  inertial  frame  of  reference.  The  integral 
form  of  the  equations  are  given  below: 
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where: 


*  V  is  the  fluid  clement  volume 

*  S  is  the  fluid  element  surface  area 

*  n  is  the  outward-pointing  unit  normal  to  S 

*  i.j.  k  are  the  Cartesian  unit  vectors 


89 


Q  = 


p 

pu 

pv 

pw 

pet 


f  = 


p) 


9  = 


h  = 


pu 
{pu2 
puv 
puw 

u(pet  +  p) 

pv 

puv 

{pv2  +  p) 
pvw 

v(pet  +  p) 

pw 

puw 

pvw 

{pw2  +  p) 
w{pet+p) 


r  ~ 


0 

rT 


s  — 


■Ty 


UTr 


'yy 

ry- 


+  vrxy  +  wrx 


UTXV  +  VTvy  +  WT, 


t  = 


0 


zy 


UTZX  +  VTzy  +  WT, 


+  AST* 


+  *T» 


90 


where: 


•  p  is  the  density 

•  p  is  the  pressure 

•  u,  v,  w  arc  the  velocity  components 

•  e  is  the  specific  energy  per  unit  volume 

•  T  is  the  temperature 

•  k  is  the  thermal  conductivity 

•  r  represents  shear  stress  components 

The  ideal  gas  laws  close  the  set  of  equations  and  the  equations  are  non-dimensionalizcd  by  freestream 
density  and  speed  of  sound.  The  turbulence  models  are  solved,  decoupled  from  the  main  conser¬ 
vation  equations  and  the  semi-discrete  form  of  the  equations  can  be  written  as: 


where  i  and  rn  denote  quantities  for  the  iih  cell  and  m,f‘  face  of  cell  i,  respectively,  and  N,  is  the 
number  of  faces  bonding  cell  i. 

The  Riemann  problem 

Ricmann  problems  and  their  solution  procedure  were  first  introduced  into  computational  fluid  dy¬ 
namics  by  Godunov  (60).  A  description  of  the  Ricmann  problem  for  a  non-stationary  flows  of  a 
perfect  gas  is  provided  here.  In  this  case,  the  state  U  of  a  perfect  gas  is  specified  completely  by 
three  dependent  variables  and  two  constants  that  are  particular  to  the  gas  under  consideration.  The 
dependent  variables  are  the  pressure,  density  and  flow  velocity,  while  the  constants  are  the  ratio 
of  the  specific  heats  and  the  gas  constant.  It  must  be  remembered  that,  all  other  state  properties 
can  be  obtained  from  the  aforesaid  variables  and  constants.  Considering  two  discrete  initial  states 
U(x,,  t} )  and  U(xi+j,  t3)  at  adjacent  nodes  at  time  level  tj  in  a  numerical  grid,  the  Ricmann  problem 
can  be  defined  as: 


U(x,t} )  =  t/j[x  <  x,  +  D(x,  —  x,+1)] 


(81) 


U{x,  tj)  -  Vi+][x  >  x i  +  Q{xi  -  x,+i)j 


(82) 
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where  JT1  is  a  number  between  0  and  1  that  fixes  the  location  of  discontinuity  between  these  two 
states.  Therefore,  the  initial  data  is  piecewise  constant,  equalling  U(.t,,  orU(x)  +  1.  t3)  depending 
on  the  location  of  the  discontinuity  between  the  two  nodes.  With  time,  the  discontinuity  between 
the  two  initial  slates  break  into  leftward  and  rightward  moving  waves,  separated  by  a  contact  sur¬ 
face  (this  depends  on  the  initial  data,  based  on  which  a  wave  may  either  be  a  shock  or  a  rarefaction). 
The  available  combinations  produce  four  unique  wave  patterns  which  are  self-similar,  depending 
only  on  x/t.  The  problem  of  determining  the  types  of  waves,  their  strengths,  and  the  flow  prop¬ 
erties  in  each  region  between  the  waves  and  contact  surface  for  some  particular  set  of  initial  data 
is  called  a  Ricmann  problem,  and  the  algorithm  for  determining  the  solution  is  called  a  Riemann 
Solver.  Since  this  is  a  one-dimensional  case,  the  pressure,  density  and  velocity  of  the  wave  on 
either  sides  of  the  contact  surface  can  be  determined  using  the  three  conservation  laws  of  mass, 
momentum  and  energy  and  the  equation  of  state.  The  conservation  laws  reduce  to  the  Rankine- 
Hugoniol  relations  for  a  shock  wave  and  an  isentropic  characteristic  equation  across  rarefaction 
waves.  These  equations  are  used  to  jump  across  the  moving  waves  into  the  unknown  state  across 
the  contact  surface.  By  assuming  that  a  single  state  exists  between  the  waves  on  the  left  and  right 
sides,  the  Ricmann  problem  may  be  reduced  to  a  single  nonlinear  algebraic  equation  in  one  un¬ 
known  for  any  particular  wave  pattern.  This  equation  is  implicit  in  pressure  or  velocity  and  can  be 
solved  iteratively. 


Rk’niann  Solver  of  Gottlieb  and  Groth 

An  efficient  Riemann  solver  will  involve  fewer  mathematical  operations  for  the  entire  solution 
procedure,  inclusive  of  the  initial  guess,  equations  used  in  the  iterative  procedure,  check  used  to 
stop  (he  iterative  procedure  and  additional  expressions  required  to  completely  specify  the  states 
on  both  sides  of  the  contact  surface  and  the  determination  of  the  wave  speeds.  From  experience, 
Gottlieb  and  Groth  found  that  numerical  computations  are  more  efficient  when  the  states  at  the  grid 
nodes  arc  defined  by  the  variables  (p,  c,  u,  7,  R)  instead  of  (p.  p,  u,  7,  R).  as  the  speed  of  sound  is  a 
dependent  variable  which  appears  more  frequently  than  the  density  in  the  equations  (60).  Instead 
of  solving  for  the  common  pressure  of  the  intermediate  states  (p*),  the  common  flow  velocity  of  the 
intermediate  states  (u*)  is  used  to  iterate  in  this  Riemann  solver,  and  the  pressure  difference  (p/*  - 
pr‘)  across  the  contact  surface  was  made  zero.  Earlier  flow  solvers  differ  from  the  Riemann  solver 
of  Gottlieb  and  Groth.  as  they  involved  iteration  with  the  variable  p\  which  was  detrimental  to 
the  computing  efficiency.  Newton’s  iterative  procedure  was  used  to  determine  the  solution  of  the 
Riemann  problem,  successive  iterates  of  u‘  are  given  by: 


uJ+i  = 


.  Pi  Ui  -  Pr  Uj 

*'  -  Pr  *'«,■*' 


(83) 


where  the  prime  denotes  differentiation  with  respect  to  u‘.  Convergence  is  determined  by  satis¬ 
fying  the  inequality  |1  —  p,*/pr’ j  <  t,  and  the  tolerance  (  is  usually  10~s.  Shock  and  rarefaction 
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equations  yield  p(*,  pr*,  C/*,  cr*,  dpi" /du'  and  dpr’/du‘.  The  relevant  equations  in  the  case  of  a 
leftward  moving  shock  wave  (u*  <  iq)  are: 


W  = 


7(  +  1  u*  —  ut 
4  ci 


(84) 


Pt‘  =  Pi  +  C,(u  -  ut)W( 


(85) 


2  CtW? 
Pl  (1+W,2) 
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ci  —  cM 


f(7t  + 1)  +  (7;  -  1  )pr_[vi 

(7i  +  1)  +  (7(  -  l)pt/pt* 


(87) 


where  Ci  —  'yipt/ct  is  computed  before  the  start  of  the  iterative  process,  and  Wi  is  the  shock 
Mach  number  with  respect  to  the  moving  gas  ahead  of  the  shock.  U)  is  a  by-product  of  the 
iterative  process.  The  shock  velocity  is  denoted  as  u;  +  C1W1  in  the  laboratory  reference  frame. 
The  variables  q*  and  V\  are  not  needed  in  the  iterative  process,  and  therefore  are  computed  later. 
Similarly,  we  could  write  a  set  of  equations  for  a  right  moving  shock  wave. For  a  left  moving 
rarefaction  wave  (u*  >  U]),  the  relevant  equations  are: 


cr  =  ct 


(u  -  Ui) 


(88) 


Pi *  =  Pi 


(89) 


pr  ~  7i  .  (90) 

Cl 

In  this  case,  q*  is  the  by-product  of  the  iterative  procedure  and  the  velocity  of  the  expansion  wave 
head  is  u/‘  +  Q*.  The  initial  guess  for  the  flow  velocity  u0*  is  obtained  by: 


Uo  — 


U/2  +  Ur 
1  +  2_ 


(91) 
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ut  =  ut  H - -Ct 

7 1  —  1 
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Ur  —  Ur  “t“  ,  Cr 

7r  —  1 


5  = 


7r  “  1  Of  Pr 
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a  =  7 liPl  >  Pr) 


(95) 


O’  =  7r(P(  <  Pr) 


(96) 


Glimm's  method 

Glimm’s  method  is  a  method  for  constructing  approximate  solutions  to  systems  of  hyperbolic  con¬ 
servation  laws  by  sampling  explicit  wave  solutions  (2).  It  is  based  on  the  idea  of  representing  a 
piecewise  continuous  flow  by  a  series  of  discontinuities.  These  waves  propagate  with  different 
speeds  and  have  different  strengths.  The  discontinuities  are  modeled  as  a  series  of  Riemann  prob¬ 
lems,  separated  by  jumps.  If  the  timcstep  is  small,  then  by  finite  propagation  speed,  the  waves 
from  adjacent  discontinuities  do  not  intersect  each  other  and  the  solutions  to  the  adjacent  Riemann 
problems  fit  together. 

Initial  conditions  for  the  Riemann  problem 

The  inviscid  fluxes  are  obtained  assuming  a  linear  data  variation  within  each  cell,  for  second  order 
accuracy.  The  left  initial  state  for  the  Riemann  problem  at  face  J  is  given  by: 

=  9i  +  f\j.sjqt  (97) 

where  q,J  -  Estimated  value 

y  r/,  -  Gradient  vector  for  cell  i 

f/  -  Vector  from  the  centroid  of  cell  i  to  the  centroid  of  cell  J  Similarly,  the  right  initial  state 
can  be  found  out.  The  gradient  vector  for  cell  i  is  found  by  a  least  squares  solution  to  the  above 
equation.  Replacing  the  unknown  value  at  the  face  centroid  with  the  known  value  at  the  centroid 
of  cell  j,  and  then  considering  all  the  nearest-neighbour  cells  of  cell  i,  we  get  the  matrix  relation: 

=  Qm~%  (98) 

In  the  above  relation,  v<(7>  *s  called  the  central  difference  gradient  as  it  is  based  upon  all  the 
nearest  neighbour  cells  of  i.  A  is  over- determined  as  a  result  of  each  cel!  having  more  nearest 
neighbours(equations)  than  unknowns.  The  solution  is  obtained  by  QR  factorization  instead  of  an 
LU  decomposition  as  QR  factorization  is  more  stable  especially  for  the  high  aspect  ratio  cells  of 
the  boundary  layer.  Constructing  the  gradient  requires  a  matrix-vector  multiplication  yielding, 

n, 

/v  *  =  -  <?,)  (99) 

m=l 
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where  w™  —  R~1Q\m  is  the  ieast-squares  weights  vector  of  cel  I  i  associated  with  nearest-neighbour 
cell  m  as  a  result  of  QR  factorization  of  A.  The  least-squares  weights  vectors,  uj™  are  computed 
only  once  and  stored  as  they  are  grid-related  quantities.  A  central  difference  gradient  is  used  for 
extremal  values  at  the  cell  centroids.  In  case  of  a  non-extremal  case,  a  one-sided  difference  least- 
squares  gradient  is  used.  A  Total  Variation  Diminishing(TVD)  scheme  is  used  to  limit  the  extremes 
at  the  cell  faces.  In  regions  of  extrema,  one-sided  differencing  is  used  resulting  in  a  “fully-upwind” 
scheme.  The  left  and  right  values  of  a  given  fluid  quantity  are  computed  for  every  face  in  the  grid, 
and  a  celt  producing  a  local  extrema  is  flagged.  A  one-sided  difference  gradient  is  then  computed 
for  each  face  of  every  flagged  cell.  The  fluid  quantity  at  each  face  oflhe  flagged  cel!  is  recomputed 
using  the  one-sided  difference  scheme  for  the  particular  face.  If  the  recomputed  fluid  quantity 
remains,  or  becomes  extremal,  it  is  limited,  by  setting  a  value  which  is  the  highest  of  the  left  and 
right  cell  centroid  values.  A  one-sided,  least-squares  gradient  can  then  be  obtained  through  a  sim¬ 
ple  correction  to  the  original  central  difference,  least-squares  gradient.  The  basic  concept  is  to 
correct  the  central-difference  gradient  in  i  by  adjusting  the  state  in  j  such  that  the  error  over  the  fit 
is  minimized.  We  have. 


where: 


A 

kj 


Cm  —  Qm  {tfi  T  t’ i  -9 >) 


( 1 00} 
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where  V  is  the  set  of  all  flagged  celts,  em  is  the  error  of  fit  at  all  the  neighbouring  cells  using  the 
central  difference,  and  r*m  is  the  vector  connecting  the  centroid  of  one  cell  with  the  centroid  of  the 
neighbouring  cell.  This  results  in: 


—  oJ 

V  9* 


>  ATi  M  } 
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where  V  q,  is  the  one-sided  least-squares  gradient  in  ceil  i  associated  with  the  face  J  and  Jj,3 
is  the  least-squares  weights  vector  of  cell  i  associated  with  nearest  neighbour  cel!  j.  If  the  near¬ 
est  cell  ncoghbours  fail  to  produce  a  change  in  the  values,  the  resulting  set  of  U  is  regarded  as 
incomplete.  This  poses  a  problem  as  the  correction  procedure  described  above  assumes  only  one 
nearest-neighbour  cell  must  be  eliminated  from  the  central  difference  stencil.  This  problem  can 
be  corrected  by  eliminating  the  contribution  of  the  flagged  cells  from  the  central  difference  fluxes 
before  applying  the  corrections.  Limiting  is  vital  as  it  ensures  stability  and  for  this,  dissipative 
effects  should  be  avoided.  This  makes  first  order  schemes  better  than  second  order  schemes.  The 
Riemann  problem  is  solved  by  a  combination  of  the  approximate  Riemann  method  of  Colei  la  (2) 
and  Newton  ’s  procedure  of  Gottlieb  and  Groth  (60). 
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Viscous  Fluxes 


The  formulation  of  the  viscous  terms  should  have  two  attributes  -  it  must  be  conservative  and  it 
must  satisfy  the  discrete  maximum  principle.  The  ratio  of  the  eddy  viscosity  to  the  molecular 
value  has  high  values  in  certain  regions  of  the  grid.  This  means  that  in  these  regions,  the  viscous 


fluxes  dominate  the  inviscid  fluxes.  Viscous  fluxes  are  required  to  add  to  the  stability  of  the  flow 


solver.  Eventhough,  they  are  expensive  and  they  account  for  a  huge  fraction  of  the  computational 
time,  their  calculation  is  worth  the  effort  due  to  the  reliability  and  convergence  properties.  The 
formulation  of  the  viscous  terms  should  be  conservative  and  must  satisfy  the  discrete  minimum 
principle.  Conservation  can  easily  be  satisfied  by  constructing  viscous  fluxes  for  the  faces,  but 
simultaneously  satisfying  the  other  conditions  is  difficult.  Efforts  are  still  on  to  construct  a  gradient 
for  the  viscous  fluxes  that  would  satisfy  all  the  conditions.  Current  studies  use  a  reasonably  good 
assumption  which  has  yielded  satisfactory  results  (36).  To  accurately  capture  the  gradients  in  the 
boundary  layers,  the  cells  are  compressed  in  the  direction  of  the  velocity  gradient. 

Temporal  Integration 

A  brief  description  of  temporal  integration  is  outlined  below.  The  governing  equations  are  recast 


as: 


+ v./r1  +  (i  -  ouvjf + v/)r  =  o 
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where  /  is  the  flux  vector  while  the  superscripts  n  and  n+1  denote  successive  time  levels  and 
0  <  d  <  1  where  0  blends  the  temporal  integration  scheme  from  fully  explicit  (fl  —  0)  to  fully 
implicit  (0  -  1).  The  discrete  forms  of  the  temporal  derivatives  are: 


(SQ  +i  =  «u(Qn+1  -  Qn)  +Qi,a(QT1  -  Qn~‘) 

'  x  i  -  At 
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The  flux  at  tn+l  is  obtained  by  noting  that  /  =  f(Q)  and  forming  a  Taylor  series  expansion  about 

r. 


(106) 


-  /S  +  AAQ  +  0(Af2) 


(107) 


where  A  =  is  the  flux  Jacobian  matrix.  This  is  split  into  two  to  account  for  the  upstream  and  the 
downstream  signal  propagation. 


AAQ  =  A+AQ+  +  A~AQ~ 


(108) 
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Substituting  this  in  the  governing  equation  for  temporal  integration,  and  rearranging  we  get  an 
equation  in  matrix  notation  as : 


LHS{Qn+ 1  -  Qn)  =  RHS  ( 1 09) 

Computation  ofthe  analytical  viscous  and  inviscid  flux  Jacobians,  due  to  the  nature  of  the  Rcimann 
solver  is  expensive.  So  the  Jacobians  of  the  Van  Leer  splitting  methods  are  used.  The  splitting  in¬ 
volves  a  temporal  damping  term  which  tends  to  make  LHS  diagonally  dominant.  The  flux  Jacobian 
matrix  A  is  split  as 

A±  =  A±rp\max  (no) 

where  [Amaj.  =  V  ■  h  +  c],  V  is  the  velocity  vector,  c  is  the  speed  of  sound,  (j  is  the  temporal 
damping  term.  The  viscous  fluxes  ensure  robustness  and  are  split  using  a  simplified  eigenvalue  ap¬ 
proach,  Temporal  damping  may  be  used  when  large  time  steps  are  used,  but  this  leads  to  a  decrease 
in  the  temporal  accuracy  and  therefore  a  sufficiently  small  timestep  is  uesd  so  that  the  temporal 
damping  is  very  less.  However,  the  addition  of  temporal  damping  to  the  viscous  fluxes  leads  to  an 
improvement  in  the  stability  of  the  solution.  This  matrix  is  solved  iteratively  using  a  symmetric 
Gauss-Seidel  procedure.  Newtons  sub-iterations  have  been  used  to  increase  the  temporal  accuracy 
of  the  unsteady  problems.  This  implicit  scheme  is  stable  for  any  Courant-Freidrichs-Lewy(CFL) 
number,  but  generally  a  value  of  106,  is  used  as  recommended  by  Strang  et  al.  (36). 

Parallelization  of  the  code 

Originally,  Cobalt  was  written  in  FORTRAN  77,  but  was  then  rewritten  in  FORTRAN  90  to  take 
advantage  of  the  new  features  such  as  runtime  usage  of  memory  and  parallelization  of  the  code. 
In  the  new  code,  it  is  also  possible  to  reduce  memory  requirements  when  the  user  requests  simpler 
physics  or  simpler  numerics.  Using  FORTRAN  90  significantly  improves  software  modularity  and 
maintainability. 

The  parallel  version  is  based  on  the  domain  decomposition  ofthe  grid,  according  to  which  each 
processor  operates  on  a  subsection  (zone)  ofthe  original  grid.  Message  Passing  Interface  (MPI)  is 
used  to  pass  information  between  the  processors.  A  variety  of  MPI  features  such  as  non-blocking 
communications,  persistent  communication  requests  and  vector  MPI  datatypes  are  used  to  mini¬ 
mize  the  communication  overhead.  The  zonal  boundaries  are  computed  exactly  the  same  manner 
as  the  zonal  interiors.  At  the  zonal  boundaries,  the  conserved  variables  for  each  boundary  cell, 
the  initial  conditions  for  the  Riemann  problem  for  each  boundary  face  and  the  gradients  in  each 
boundary  cell  (viscous  cases  only)  are  passed  to  the  processors. 

A  preprocessing  program  is  used  to  partition  the  grid  into  zones.  The  intermediate  decomposi¬ 
tion  is  done  by  METIS  (15).  This  helps  in  creating  almost  equally  balanced  zones  in  such  a  way 
that  the  sizes  of  the  zones  are  almost  identical  and  excellent  load  balancing  is  achieved  between 
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the  processing  nodes,  A  block  Gauss-Siedcl  scheme  is  used  along  with  sub-iterations  to  increase 
the  temporal  accuracy  of  unsteady  problems. 
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